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INTRODUCTION 


The  major  objective  of  this  report  is  to  establish  a  relationship 
between  fatigue  crack  growth  and  internal  damping.  In  defining  this 
relationship,  a  basis  for  monitoring  the  structural  integrity  of  systems 
subject  to  fatigue  will  be  demonstrated.  Thus,  the  value  of  this  rela¬ 
tionship  is  in  nondestructive  diagnostics,  where  the  detection  of  crack 
initiation  and  propagation  are  of  great  interest. 

A  secondary  objective  is  the  mathematical  modeling  of  the  crack 
growth-internal  damping  relationship.  To  understand  the  relationship, 
and  apply  it  in  the  analysis  of  complex  structures,  a  model  is  essential. 
Without  a  model,  extension  of  the  relationship  to  allow  monitoring  of 
structural  integrity  would  be  awkward  and  prone  to  failure.  The  model 
provides  a  foundation  to  begin  observations  of  fatigue  failure,  in 
complex  structures,  using  damping. 

Previous  work  that  categorized  and  analyzed  damping  mechanisms  and 
explored  their  relation  with  fatigue  is  described.  Also,  a  detailed 
description  of  subgoals  required  to  establish  a  crack  growth-internal 
damping  relationship  is  provided. 


BACKGROUND 

The  purpose  of  this  background  section  is  to  review  existing 
literature  about  internal  damping  and  crack  propagation  and  give  the 
reader  an  introduction  to  the  research  topic.  Aspects  of  damping  or 
fatigue  that  relate  to  the  research  objective  and  are  crucial  to  its 
accomplishment  are  explored. 


Damping  Mechanisms 


When  a  solid  is  subjected  to  vibration,  several  processes  may  occur 
to  dissipate  energy,  even  at  low  amplitudes  of  vibration.  The  ratio  of 
energy  dissipated  to  peak  elastic  energy  is  usually  quite  small;  the 
main  significance  being  that  information  is  provided  as  to  the  nature  of 
the  processes  involved.  At  large  amplitudes,  the  proportion  of  energy 
dissipated  may  be  much  larger  as  permanent  damage  accumulates  within  the 
material . 

The  most  important  material  damping  process  is  the  vibration  of 
dislocations.  A  dislocation  is  defined  as  a  stable  arrangement  of  atoms 
such  that,  in  the  region  of  a  few  atomic  distances,  (n  +  1)  atoms  in  the 
slip  plane  face  (n)  atoms  across  the  slip  plane  (Figure  1).  Movement  of 
dislocations  produce  an  incremental  strain  that  diminishes  the  effective 
elastic  modulus,  while  the  energy  dislocations  dissipate  during  vibration 
manifesting  itself  as  "internal  friction."  There  are  also  several  relaxa¬ 
tion  processes  contributing  to  internal  friction  that  do  not  involve 
dislocations.  They  are  called  relaxation  processes  because  they  change 
the  shape  of  the  solid  under  load  and  because  each  has  a  certain  charac¬ 
teristic  "relaxation  time."  Such  processes  are  grain  boundary  sliding, 
thermal  diffusion,  atomic  diffusion,  and  magnetostriction. 

The  relaxation  processes  were  named  anelastic  by  Zener,  who  first 
interpreted  them  in  toto  (Ref  1).  They  create  a  peak  in  the  curve  of 
internal  friction  versus  frequency  of  vibration  at  a  frequency  approxi¬ 
mately  equal  to  (relaxation  time)  .  The  magnitude  of  the  internal  fric¬ 
tion  is  independent  of  the  amplitude  of  vibration  in  a  first  approximation. 
The  more  significant  dislocation  process  is  different  in  that  no  peak 
exists  until  megacycle  frequencies  are  reached,  and  the  magnitudes  of 
the  associated  internal  friction  are  not  so  independent  of  the  amplitude 
of  vibration.  To  differentiate  it  from  the  relaxation  processes  it  has 
been  called  a  process  of  static  hysteresis  (Ref  2). 

Internal  friction  typically  varies  with  the  amplitude  of  vibration 
(Ref  3)  (Figure  2).  Up  to  a  strain  amplitude  of  about  10  ^  (Figure  2) 
internal  friction  is  independent  of  amplitude  ("amplitude-independent" 
range),  then  commences  to  increase  with  amplitude  ("amplitude-dependent" 


range).  For  very  large  amplitudes,  material  enters  the  fatigue  range 
and  damping  increases  still  faster,  ultimately  resulting  in  a  fatigue 
fracture.  In  addition  to  the  amplitude  of  vibration,  temperature  and 
frequency  may  affect  the  internal  friction. 

Considerable  evidence  indicates  that  most  of  this  internal  friction 
is  connected  with  the  movement  of  dislocations  (Ref  4  through  8).  Pro¬ 
vided  the  selected  conditions  do  not  excite  a  relaxation  peak,  dislocation 
movement  is  the  process  responsible  for  most  of  the  internal  friction  in 
metals . 

Having  established  the  basics  of  internal  damping  and  its  primary 
mechanisms,  the  next  step  is  to  examine  previous  attempts  to  relate 
damping  and  fatigue. 

Damp ing  and  Fatigue 

Various  researchers  have  sought  a  relation  between  damping  and 
fatigue.  For  many  years,  damping  was  believed  to  be  the  first  indica¬ 
tion  of  damage  to  material  subjected  to  cyclic  stress.  It  was  further 
believed  that  each  metal  was  capable  of  dissipating  only  a  certain 
quantity  of  energy  by  internal  damping,  and  that  when  this  quantity  was 
exceeded,  fracture  by  fatigue  must  follow.  In  fact,  every  experiment 
made  for  the  purpose  of  verifying  this  hypothesis  has  given  negative 
results . 

Experimentally,  Foppl  (Ref  9)  has  shown  that  the  capacity  for 
dissipation  of  energy  in  the  form  of  heat  by  internal  damping  is  un¬ 
limited.  He  conducted  an  alternating  torsion  test  continuously  for  a 
period  of  3  years  without  the  test  piece  breaking  or  showing  any  sign  of 
deterioration  after  1,100  million  cycles  of  stress.  The  rate  of  dissi¬ 
pation  of  energy  was  sufficient  to  maintain  the  test  piece  at  about 
100°C  above  the  ambient  temperature. 

Damping,  then,  is  not  directly  related  with  fracture  by  fatigue. 

It  depends  principally  on  the  strain  energy  of  the  cycles,  on  the 
temperature  and  certain  other  conditions  of  the  test,  and  also,  to  a 
smaller  extent,  on  the  number  of  stress  cycles  which  have  previously 
been  applied.  For  many  metals,  the  damping  capacity  increases 
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noticeably  up  to  a  certain  number  of  cycles  (about  2  to  10  million),  and 
then  remains  stationary  as  long  as  the  conditions  of  the  endurance  test 
are  maintained  unchanged. 

From  previous  tests,  it  is  possible  to  advance  the  following  theory: 
a  fatigue  crack  results  from  the  accumulation  of  damping  effects  (i.e., 
from  the  fact  that  under  a  variable  stress  higher  than  the  fatigue  limit, 
constantly  applied,  damping  continually  increases  up  to  fracture).  For 
a  cyclic  stress  below  the  fatigue  limit,  the  damping  moves  toward  a  fixed 
value,  or  else  diminishes  continuously. 

The  origination  of  fatigue  cracks  is  marked  by  an  increase  in  damping. 
This  increase  is  associated  with  the  accumulation  of  dislocations.  Since 
cyclic  input  of  strain  energy  alters  the  thermal  and  internal  surface 
energies,  dislocations  tend  to  collect  on  the  planes  already  containing 
one  or  more  dislocations  (Figure  3).  Eventually,  (n)  atoms  face  (n  +  2), 

(n  +  3),  or  more  atoms,  resulting  in  relative  motion  of  one  plane  to  the 
other.  For  aluminum  alloys,  Hanstock  and  Murray  (Ref  10)  were  able  to 
predict  the  endurance,  in  tests  under  alternating  torsion,  by  observing 
the  stage  at  which  this  increase  in  damping  occurred. 

The  change  in  damping  with  the  number  of  stress  cycles  applied  is 
related  to  the  process  of  accumulative  damage  that  results  in  fatigue 
fracture,  although  it  does  not  directly  express  the  effect  of  damage, 
but  rather  the  effect  of  localized  inelastic  deformations  accompanying 
damage.  It  is,  therefore,  not  the  absolute  value  of  damping  capacity, 
but  its  change  with  the  number  of  cycles  that  indicates  the  process  of 
progressive  damage.  As  a  further  complication  to  the  interpretation  of 
damping,  above  a  certain  value  of  stress,  the  damping  capacity  is  reduced 
by  work  hardening,  and  is  increased  by  the  repetition  of  stress  cycles. 

The  work-hardening  process  in  a  material  seems  to  be  due  to  the 
long-range  elastic  Interactions  of  its  dislocations.  Additional  input 
of  strain  energy  may  help  to  overcome  these  long-range  stresses  to 
accumulate  more  dislocations  on  the  slip  planes,  until  the  crystal 
separates  into  two  or  more  fragments,  with  interference  and  distortion 
of  the  adjacent  crystals.  Continued  cycling  would  eventually  cause  the 
cracks  to  grow  and  coalesce  into  an  open  area  (resulting  in  an  insuffi¬ 
cient  area  to  support  the  load)  and  finally  fracture. 
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With  the  lack  of  success  in  previous  attempts  to  directly  relate 
crack  growth  and  internal  damping,  a  different  approach  was  sought. 
Prerequisites  to  defining  a  new  approach  are  the  examination  of  fatigue 
processes  and  an  understanding  of  how  these  processes  are  modeled. 

Fatigue  Processes 

In  examining  the  crack  tip  on  a  microscopic  scale,  it  is  convenient 
to  regard  the  fatigue  life  as  consisting  of  a  number  of  stages  (Ref  11, 
12).  The  total  number  of  cycles  to  fracture  is  divided  into  four  parts. 

ic 

The  first  stage  occupies  an  average  period  of  cycles  and  represents 

the  completion  of  work  hardening.  A  plastic  zone  forms  under  load  and 

* 

the  crack  tip  blunts.  An  average  period  N_  is  required  for  the  formation 

of  the  first  submicrocracks,  created  by  the  accumulation  of  dislocations 

at  slip  planes.  Further  cycles  delineate  the  third  stage,  during 

which  the  submicrocracks  grow  and  coalesce  to  form  a  crack  of  detectable 

size.  Finally,  these  cracks  propagate  across  grains  until  fracture  or 

* 

rupture  occurs,  and  an  additional  number  of  cycles  elapse  denoting 
the  fourth  stage. 

Damage  occurs  in  the  second  stage  and  can  be  attributed  to  sub¬ 
microcracks  that  are  formed  in  regions  where  stresses  are  concentrated 
because  of  the  piling  up  of  dislocations,  or  the  presence  of  inclusions. 
In  succeeding  stages  two  possible  developments  can  occur.  Either  the 
cracks  grow  quickly  to  a  detectable  size  (the  observable  cracks  at  the 
end  of  the  third  stage)  or  plastic  flow  occurs  at  the  crack  ends  which 
are  stress  concentration  sites,  so  that  the  crack  growth  is  temporarily 
stayed.  Stresses  subsequently  increase  because  of  work  hardening,  and 
new  submicrocracks  are  formed  at  the  ends  of  old  submicrocracks.  This 
would  suggest  that  cracks  grow  iteratively. 

The  general  and  most  significant  features  of  fatigue  failure  are 
the  initiation  of  surface  microcracks  and  their  subsequent  extension 
across  and  penetration  into  the  body  of  the  metal  The  increased  life 
resulting  from  the  removal  of  a  surface  layer  at  frequent  intervals 
throughout  a  test,  irrespective  of  whether  the  life  is  many  millions  or 
only  a  few  thousand  cycles  (Ref  13),  demonstrates  that  crack  initiation 
is  confined  to  the  surface  grains.  Dislocation  models  leading  to  surface 
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cracking  have  been  proposed  (Ref  14,  15,  and  16)  that  lead  to  a  geometric 
cause  of  damage  or  damage  on  the  slip  plane.  In  general,  surface  cracking 
can  occur  in  any  of  three  ways;  specifically,  as  a  continuation  of  surface 
roughening  in  broad  slip  bands,  as  a  result  of  severe  strain  incompati¬ 
bilities  across  grain  boundaries,  or  because  of  the  presence  of  inclusions 
or  inhomogeneities  in  the  surface. 

Theories  of  crack  growth  based  on  a  dislocation  model  begin  with  a 
crack  loaded  in  antiplane  strain  (Mode  III).  The  plastic  zone  at  its 
tip  can  be  conveniently  represented  by  a  continuously  distributed  array 
of  infinitesimal  dislocations  on  the  crack  plane.  Crack  growth  is  assumed 
to  start  when  the  accumulated  plastic  strain  distribution  at  a  crack  tip 
exceeds  a  critical  value  (stage  3),  and  continues  as  this  value  is  exceeded 
at  successive  points  ahead  of  the  original  crack  tip.  It  is  then  assumed 
that  the  behavior  in  Mode  I  is  similar  to  behavior  in  Mode  III.  In  general, 

such  theories  (Ref  17)  predict  that  the  rate  of  crack  growth  is  proportional 
4 

to  Kj  (or  AFC,  depending  on  the  load  cycle  considered,  where  AK  =  (KmflX 
-  K  .  )  and  K  and  K  .  are  the  maximum  and  minimum  values  of  KT  during 
the  fatigue  cycle).  As  an  example,  Weertman's  theory  (Ref  18)  started 
with  a  model  (Ref  19)  of  a  freely  slipping  crack  subjected  to  an  applied 
shear.  In  the  model,  a  crack  of  length  2a,  having  plastic  zones  at  either 
end,  was  considered  and  the  dislocation  theory  was  used  to  calculate  the 
displacement  in  the  vicinity  of  the  crack  tips  due  to  this  shear  stress. 

It  was  then  shown  that,  under  certain  conditions,  a  crack  lying  normal 
to  a  tensile  stress,  o,  can  be  considered  in  a  similar  manner. 

The  energy  associated  with  the  plastic  zone  at  a  crack  tip  is  pro- 

4 

portional  to  .  Theories  based  on  the  energy  required  to  operate  the 

fracture  mechanism,  in  general,  will  predict  that  the  rate  of  crack  growth 

4 

is  proportional  to  ,  which  is  in  agreement  with  dislocation  theories. 

The  modeling  of  fatigue  processes  has  found  some  success  in  estimating 
fracture  based  on  energy  or  work  at  the  crack  tip.  In  particular,  equa-  > 

tions  relating  plastic  work  and  fracture  are  comi  on  in  the  literature. 

Relating  crack  growth  to  internal  damping  through  plastic  work  dictates 
a  study  of  the  role  of  plastic  work  in  fatigue. 
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Plastic  Work  in  Fatigue 


Coffin  and  others  (Ref  20)  found  that  data  obtained  from  reversed 
direct  constant  plastic  strain  amplitude  tests,  designed  to  6ive  failures 
in  up  to  about  10^  cycles,  conformed  to  the  relationship: 

N1/2  e  =  P  (1) 

P 

where  N  is  the  number  of  fatigue  cycles,  e^  is  the  plastic  strain  ampli¬ 
tude,  and  C  is  a  material  constant.  In  the  case  of  annealed  low-carbon 
steels  exhibiting  a  marked  yield  point,  the  relation  does  not  apply  in 
the  immediate  region  of  discontinuous  yielding,  but  becomes  applicable 
for  the  balance  of  the  stress -strain  curve.  Manson  (Ref  21  through  24) 
presented  an  analysis  of  direct  loading  constant  plastic  strain  amplitude 

tests  for  30  materials  of  widely  differing  static  properties  (Ref  25) 

4 

showing  that  in  all  cases  the  endurance  up  to  about  5  x  10  cycles  was 
related  to  the  plastic  strain  range  by  the  relationship. 


He  suggested  that  the  results  conformed  to  this  relationship,  and  not 
Equation  1,  because  N  had  been  taken  as  the  number  of  cycles  to  complete 
the  fracture  of  the  specimen.  Equation  1  was  based  on  endurances  where 
cracks  were  first  visible  on  the  specimen  surface.  It  is  now  accepted 
that  the  value  of  the  exponent  in  Equation  1  is  not  a  universal  constant; 
it  varies  with  material,  environment,  and  the  criterion  adopted  for  failure 
(Ref  26,  27). 

The  relationship  between  stress,  plastic  strain,  and  plastic  strain 
energy  during  a  low-endurance  fatigue  test  was  studied  by  Morrow  (Ref 
26).  He  related  fatigue  life  to  the  total  plastic  strain  energy  generated; 
that  is,  C  is  made  a  function  of  the  integrated  areas  of  the  hysteresis 
loop  up  to  the  point  of  failure.  Radhakishnan  (Ref  28)  modeled  fatigm 
failure  on  the  basis  of  total  energy  absorbed  and  proposed  that  the  total 
hysteresis  energy  to  failure,  W^.,  can  be  given  by: 
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wf  ■  / 


dw(N)  dN 


where  dw  is  the  plastic  energy  absorbed  per  cycle  at  any  given  time  and 
Nj.  is  the  number  of  fatigue  cycles  to  failure.  For  most  materials,  the 
hysteresis  loop  is  relatively  stable  after  a  certain  initial  life.  Taking 
e  as  the  average  value  of  the  plastic-strain  range  during  the  lifetime, 
the  energy  absorbed  per  cycle  can  be  approximately  given  as: 


dw  =  A„o  e 
2  a  p 


and  the  total  energy  absorbed,  W^,  up  to  fracture  will  be: 


W.  =  A.o  e  N 
f  2  a  p  f 


where  A  is  a  mater i al -dependent  constant  and  0  is  the  maximum  cyclic 
stress . 

4 

Experimental  data  supporting  the  (AK)  depe  uency  of  growth  rate 
predicted  by  both  dislocation  and  energy  theories  can  be  found  in  the 
literature.  For  example,  tests  (Ref  29)  on  steel,  aluminum  alloy, 
titanium,  and  magnesium  sheet  have  shown  that: 


D(  AK) 


where  I)  is  a  material  constant.  Similar  results  were  obtained  for  70/30 
brass  specimens  (Ref  30).  On  the  other  hand,  data  can  be  found  in  the 
literature  that  do  not  conform  to  this  relationship  (Ref  31). 

Attempts  to  correlate  plastic  work  directly  with  crack  growth  (Ref 
32)  have  resulted  in: 
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where  A  is  a  dimensionless  constant,  G  is  the  shear  modulus,  if>  is  a 
strength  parameter,  and  U  is  the  integrated  plastic  work  expended  in  the 
plastic  zone.  Izumi  and  Fine  (Ref  32)  conclude  that  the  plastic  work 
per  unit  area  of  fatigue  crack  propagation  is  an  important  parameter 
that  establishes  the  rate  of  fatigue  crack  growth.  This  report  will 
assume  that  Izumi  and  Fine  are  accurate  in  their  conclusion,  and  will 
attempt  to  use  plastic  work  as  a  basis  for  measuring  structural  damage 
and  coupling  internal  damping  with  crack  propagation. 


OBJECTIVES 

This  research  has  two  objectives.  The  first  objective  is  to 
experimentally  develop  a  relationship  between  lew  cycle  fatigue  crack 
growth  and  internal  damping.  The  second  objective  is  the  mathematical 
modeling  of  specimens  to  interpret  and  extrapolate  the  internal  damping- 
crack  growth  relationship  and  apply  it  to  other  geometries  and  load 
condit ions . 

Accomplishment  of  these  two  objectives  will  establish  a  foundation 
upon  which  nondestructive  evaluation  techniques  may  be  constructed.  The 
first  objective,  the  relationship,  provides  a  basis  for  monitoring 
structural  integrity.  The  second  objective,  the  mathematical  model, 
provides  a  process  for  expanding  the  relationship  to  meet  the  demands  of 
monitoring  complex  structures. 

Relationship  of  Crack  Growth  and  Damping 

Various  researchers  have  shown  that  crack  propagation  and  internal 
damping  are  closely  linked:  fatigue  damage  is  accompanied  by  changes  in 
specific  damping  capacity.  This  relationship,  however,  is  not  a  direct 
one.  The  failure  of  many  experiments  to  prove  a  direct  relationship  has 
made  that  point  clear.  To  determine  how  damping  is  related  to  crack 
growth,  an  investigation  must  look  for  an  indirect  linkage  between  the 
two  phenomenon.  An  indirect  linkage  suggests  that  the  establishment  of 
this  relationship  requires  the  merger  of  two  previously  developed  theories 
rather  than  the  creation  of  a  new  theory;  two  theories,  an  internal  damping 
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mechanism  and  a  fatigue  damage  mechanism,  that  have  a  common  element. 

In  merging  the  two  mechanisms  to  form  a  single  theory  of  fatigue  damage 
accompanied  by  internal  damping,  care  must  be  taken  to  assure  that  both 
mechanisms  are  used  within  the  context  of  their  assumptions. 

The  relationship  must  be  expressed  in  a  manner  that  facilitates 
experimental  verification.  Thus,  an  analytical  expression  of  crack 
growth  or  fatigue  damage  as  a  function  of  changes  in  specific  damping 
capacity  is  required.  To  verify  this  expression  and  the  associated 
relationship,  a  series  of  fatigue  tests  must  be  planned.  In  general, 
testing  would  begin  with  fatigue  specimens  that  minimize  effects  of 
complex  parameters  such  as  specimen  geometry,  variable  loading  ind 
stress  fields.  After  these  initial  experiments  provided  satisfactory 
results,  more  complex  parameters  could  be  included  for  testing  and 
verification. 

Mathematical  Model 

To  incorporate  increasing  complexity  into  the  testing,  the 
specimens  must  be  mathematically  modeled.  The  modeling  provides  two 
significant  functions: 

1.  A  model  allows  more  complex  specimens  to  be  fatigue  tested  by 
revealing  the  interaction  of  the  additional  parameters.  Changes  to 
specimen  geometry  or  loading  will  affect  both  crack  growth  and  internal 
damping,  the  effects  of  these  changes  must  be  known  in  advance  in  order 
to  verify  the  test  results. 

2.  A  model  will  provide  a  solid  foundation  for  a  diagnostic  tool 
for  detecting  crack  growth  in  structures.  The  model  must  be  simple  to 
act  as  a  tool,  but  it  must  be  comprehensive  to  provide  accuracy. 

Verification  of  the  model  requires  simple  fatigue  specimens  that 
can  be  vibrated  at  low  amplitudes  during  the  fatigue  process  for  data 
collection.  The  data  from  the  series  of  tests  will  be  analyzed  using 
the  model  for  comparison  with  the  actual  test  results. 
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The  model  will  have  to  include  elements  to  express  a  continuous 
structure  (the  specimen)  undergoing  both  dynamic  (data  collection 
cycles)  and  static  (fatigue  cycle)  loading.  In  addition,  the  model  must 
adequately  handle  a  growing  crack  (or  cracks)  in  the  specimen  as  well  as 
the  changing  internal  damping  capacity.  The  internal  damping  is  not 
only  changing  at  a  global  level,  but  on  a  local  level  as  well.  Damping 
is  a  localized  phenomenon,  and  areas  close  to  a  propagating  crack  will 
dissipate  heat  at  higher  rates  than  areas  remote  from  the  crack  (Ref 
33) . 

TECHNICAL  APPROACH 

An  analytical  approach  to  develop  a  relationship  between  crack 
growth  and  internal  damping  is  presented  to  provide  the  theoretical 
foundation  for  the  relationship;  a  foundation  integrating  internal 
damping,  crack  growth,  plastic  strain,  and  dislocations. 

Internal  Damping  Mechanism 

The  goals  of  this  section  are  to  define  an  analytical  expression 
for  internal  damping  and  to  accommodate,  within  that  expression,  the 
strain  distribution  experienced  by  a  structure  or  system.  The  phrase 
"analytical  expression"  is  used,  rather  than  "mathematical  model," 
because  the  product  of  this  section  represents  a  single  mathematical 
term  or  expression  of  the  final  model. 

Analytical  Expression.  In  1956,  Granato  and  Liicke  (Ref  34)  pre¬ 
sented  a  theory  of  internal  damping  due  to  dislocations.  Granato  and 
Lucke  examined  two  types  of  losses:  frequency-dependent  loss  (A^)  and 
strain-amplitude  dependent  hysteresis  loss  (A,,).  Their  theory  provides 
a  quantitative  interpretation  of  these  losses.  The  two  losses  are  given 
by  (Ref  35): 
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where  A  =  mechanical  energy  lost  per  cycle  divided  by  twice 
the  total  vibrational  energy 

A  =  total  length  of  movable  dislocation  line  per  unit  volume 
(dislocation  density) 

a  =  lattice  parameter 

Eg  =  strain  amplitude 

u  =  circular  frequency 

v  =  Poisson's  ratio 

AQ  =  4(1-v)/tt2 

A  =  irpb2 

C  =  2Gb2/u( 1-v) 

p  =  density 

b  =  Burger's  vector 

=  network  loop  length 

=  impurities  loop  length 

Le  =  effective  length  of  dislocation  loop 

K  =  orientation  and  anisotropy  parameter 

£2  =  orientation  factor 


Assuming  that  all  parameters  in  Equation  5  can  be  held  constant  for 
a  given  material  under  cyclic  stress,  with  the  exception  of  A  and 
the  equation  may  be  rewritten  as: 


12 


(6) 


AH  =  AC1C2/E0  exP('C2/'e0') 

where  =  RA^L^  /itL^  and  =  Ke'a/L^.  As  modeled  by  Equation  6,  all 
of  the  strain  amplitude  data  may  be  described  in  terms  of  a  single  func¬ 
tion  and  when  plotted  in  the  form: 

log(e0AH)  versus  \/tQ 

should  lie  on  a  straight  line.  The  slope  of  this  straight  line  is  given 
by  and  the  intercept  by  AC^C^-  Figure  4  shows  a  Granato-Lucke  plot 
of  strain  amplitude  times  decrement  versus  (strain  amplitude)  ^ . 

Bauer  and  Gordon  (Ref  36)  expressed  the  ampl itude- independent 
decrement,  A^.,  resulting  from  dislocation  vibration  in  the  kilocycle 
range  of  frequencies  as: 

A  i  =  C3A  (7) 

4  2 

where  =  RA^L^  Bw/h  C.  The  total  decrement,  A,  equal  to  A^  +  A^  is 
now  presented  in  terms  of  Equation  6  and  7. 

A  =  A[C3  +  C1C2/eq  exp(C2/E0)]  (8) 

Equation  8  provides  a  relation  between  internal  damping,  A,  and  dislocation 
density,  A. 

The  Granato-Lucke  (G-L)  theory  does  not  always  reflect  experimental 
results  and  its  accuracy  is  unclear.  At  times,  experimental  data  plots 
of  log(e  A  )  versus  l/en  are  not  straight  but  somewhat  curved  (see 
Figure  4).  The  decrement  measured  at  low  strain  amplitudes  appears  to 
be  greater  than  that  predicted  by  the  theory.  Some  have  argued  that  the 
model  is  essentially  correct  but  the  calculations  are  not  yet  sufficiently 
refined.  Other  authors  have  recommended  modifications  to  correct  the 
model  (Ref  37  through  40).  Unfortunately,  all  these  detailed  discussions 
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are  inconclusive  because  of  the  difficulty  of  specifying  all  the  various 
significant  factors.  In  particular,  the  dependence  of  the  decrement  on 
strain  amplitude  (Equation  5)  is  largely  determined  by  the  distribution 
function  for  the  lengths  of  the  dislocation  loops,  L^.  Granato  and  Liicke 
assume  an  exponential  function,  but  there  is  little  substantiation  of 
this  approximation  for  any  given  material  or  geometry. 

The  theoretical  treatment  that  researchers  have  used  to  explain 
internal  damping,  as  illustrated  in  Figure  2,  is  broken  into  three 
regions.  These  regions  are  labeled  in  Figure  2.  The  first  two  regions, 
elastic  and  anelastic,  are  a  function  of  dislocation  damping.  For  the 
elastic  region,  the  dislocations  are  assumed  pinned  by  point  defects, 
and  vibrating  as  stretched  strings,  subject  to  viscous  damping.  In  the 
anelastic  region,  the  stress  level  is  sufficient  to  bow  out  the  disloca¬ 
tion  and  pull  it  away  from  pinning  point  defects.  Lowering  the  stress 
repins  the  dislocation,  but  the  motion  is  different  on  the  outward  and 
return  paths,  producing  a  hysteresis  loop.  The  third  region,  fatigue- 
plastic,  is  believed  to  result  from  the  formation  of  microcracks  that 
can  fuse  together  to  form  fatigue  cracks.  For  pure  aluminum,  shown  in 

Figure  2,  dislocation  damping  is  the  dominate  form  of  damping  up  to  2  x 
-4 

10  microstrain. 

A  second  treatment,  evolving  from  McKavangh  and  Stacey  (Ref  41)  and 
Jackson  and  Anderson  (Ref  42),  is  that  the  dominate  cause  of  energy 
dissipation,  in  the  anelastic  region,  is  grain  boundary  relaxation  or 
slippage.  It  is  difficult  to  apply  the  conclusions  of  Jackson  and 
Anderson  to  this  research,  however,  since  their  study  investigated 
attenuation  in  the  earth's  mantle,  an  environment  of  high  temperature 
where  dislocations  are  likely  to  be  annealed  out.  The  work  of  McKavangh 
and  Stacey,  on  the  other  hand,  raises  doubt  as  to  the  nature  of  the 
hysteresis  in  the  anelastic  region.  Their  observations  support  non¬ 
linear  internal  friction  due  to  cracks  or  grain  boundaries.  Further 
examination  into  the  nature  of  anelastic  internal  damping  is  needed,  but 
given  the  two  theoretical  treatments  available,  this  research  effort 
will  rely  upon  the  extensive  experimentation  and  analysis  that  are 
embodied  in  the  Granato  and  Liicke  model. 
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Application  of  the  G-L  theory  in  the  form  of  Equation  8  is  compli¬ 
cated  by  the  pattern  of  strain  amplitude  distribution  within  a  specimen. 
Longitudinal,  flexural,  and  torsional  forced  vibrations  produce  different 
strain  distributions.  The  different  strain  distributions  necessitate 
the  derivation  of  mathematical  expressions  of  the  decrement  on  the  basis 
of  each  type  of  vibration. 

Strain  Distribution.  According  to  Equation  4  of  Reference  43,  the 
strain  amplitude-dependent  decrement  distribution  is  given  by: 

Ah  =  f  e2AH(E)dV/f  e2dV  (9) 

where  V  is  the  volume  of  the  specimen.  Povolo  (Ref  44)  showed  that  an 
expression  for  A.,  can  be  obtained  directly  in  terms  of  exponential 
integral  functions  and  that  Equation  6  oversimplifies  the  results  for 
the  valid  range  of  strains  as  given  by  Granato  and  Liicke,  that  is, 

E0^2  ~  0-25.  From  Povolo: 

a|J  =  ac1(c2/e0)2e1(c2/e0)  (10) 

where  E^  is  the  exponential  integral  of  order  n  =  1  and: 


En(x) 


is  the  exponential  integral  function.  From  Equation  9,  using  Povolo' s 
corrected  strain  amplitude-dependent  decrement  of  Equation  10,  the  dis- 
trihuMon  peculiar  to  the  type  of  vibration  can  be  derived. 

A  slight  modification  to  Equation  9  (required  for  specimens  that 
have  developed  interior  plastic  zones  at  a  location  of  stress  concen¬ 
tration)  is  separate  integrals  for  the  plastic  and  elastic  regions. 
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Given  a  small  amplitude  of  vibration,  e^,  it  is  assumed  that  the  only 
parameter  that  changes  in  Equation  9  between  the  two  regions  is  the 
dislocation  density,  A.  This  assumption  is  based  on  small  amplitude 
vibration  of  the  specimen  when  internal  damping  data  is  collected. 
Burdett  and  Queen  (Ref  45)  noted  that  internal  damping  developed  during 
Pi  nstic  strain  may  be  expressed  in  the  form: 

Ap  =  A  exp  [ B ( e Q  -  £p)] 


where  e  is  the  critical  strain  above  which  plastic  strain  occurs,  and  A 
P 

and  B  are  temperature-dependent  terms.  Burdett  and  Queen  concluded  that 
the  G-L  model  is  valid  in  cases  where  the  movement  of  the  dislocations 
is  not  affected  by  obstacles  outside  their  equilibrium  position.  This 
is  the  case  when  the  specimen  experiences  low  strain  vibration;  therefore, 
it  is  valid  to  assume  that  only  the  dislocation  density  has  changed  when 
comparing  internal  damping  between  the  elastic  region  and  plastic  zone. 

The  distinction  between  the  A^  and  Ap  contributions  appears  to  involve 
the  degree  of  pinning  of  dislocations.  For  lightly  pinned  dislocations 
the  contributions  are  well  separated.  In  the  case  of  heavily  pinned 
dislocations,  the  A  is  suppressed  completely  and  the  dislocations 
contribute  immediately  to  Ap  at  relatively  high  strain  amplitudes. 

The  plastic  zone  size  and  shape  depends  on  the  plastic  flow 

2 

properties  (Ref  46),  but  the  dimensions  are  proportional  to  (K^/o^)  . 

The  nominal  plastic  zone  radius,  r  ,  is  approximately: 


r 

P 


for  plane  stress 


(11) 


and 


r 

P 


1_ 

6tt 


for  plane  strain 


(12) 
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where  is  yield  stress.  The  plastic  zone  radius  is  about  half  the 
extent  of  the  plastic  zone  and  is  applicable  when  the  zone  size  is  small 
compared  to  the  specimen  dimensions  in  the  plane  of  the  plate. 

Breaking  Equation  9  into  two  separate  integrals  (the  first  for  the 
elastic  region  and  the  second  for  the  plastic  zone)  produces: 

K  =  fAn  L  E20jkE)dV  +  A  f  E20^(E)dV]/.f  E2dV  (13) 

rl  €»  e  n  p  p  H 

where  the  parameters  with  subscripts  e  and  p  correspond  to  the  elastic 
region  and  plastic  zone,  respectively;  the  strain-dependent  function, 

0^,  is  given  by: 

=  wvVw  <i4) 

C  C 

that  is,  0  =  A  /A.  In  a  similar  manner,  the  strain-independent  decrement 

H  H 

— C 

may  be  expressed  as  A.  =  C,  (A  V  +  A  V  )/V,  where  V  and  V  are  the 
r  I  3  e  e  p  p  e  p 

volumes  of  the  elastic  region  and  plastic  zone,  respectively. 

Equation  13  expresses  material  dislocation  damping  as  a  function  of 
elastic  and  plastic  strain  distributions  given  by  Equation  14.  Having 
obtained  an  analytical  expression  for  dislocation  damping,  the  second 
step  is  to  develop  a  fracture  failure  criterion. 

Crack  Growth 

Monitoring  structural  integrity  requires  a  procedure  for  quantifying 
fatigue  damage  and  estimating  crack  growth.  A  fracture  criterion  based 
on  plastic  work  is  proposed  for  accomplishing  this  task  and  providing  an 
interface  to  amplitude-dependent  damping. 

Failure  Criterion.  In  selecting  Equation  3  as  a  fracture  criterion, 
Feltner  and  Morrow  used  the  empirical  equation  0fl  =  k^e^11  (see  Figure  5) 
to  relate  the  true  plastic  strain  to  true  stress  as  shown  in  Figure  6 
(Ref  47).  With  their  empirical  stress -strain  relation,  Feltner  and 
Morrow  were  able  to  express  Equation  3  as: 
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(15) 


so  that  the  total  energy  absorbed  to  fracture  is  a  function  of  the  true 
plastic  strain  average,  e^.  For  the  purposes  of  this  investigation,  the 
definitions  of  true  stress  and  true  strain  are  taken  from  Davis,  Troxell, 
and  Wiskocil  (Ref  48).  In  postulating  that  the  plastic-strain  energy 
portion  of  the  hysteresis  energy  accounts  for  the  damaging  effects  of 
cyclic  stress,  Feltner  and  Morrow  made  three  assumptions:  (1)  a  loga¬ 
rithmic  plot  of  static  true  stress  versus  true  plastic  strain  is  valid 
when  extrapolated  back  into  the  fatigue  stress  region;  (2)  the  damaging 
energy  per  cycle  for  a  given  stress  amplitude  is  constant  and  is  equal 
to  the  area  under  the  static  stress  plastic-strain  curve;  and  (3)  the 
total  damaging  energy  required  to  cause  fatigue  fracture  is  constant  and 
as  a  first  approximation  is  equal  to  the  area  under  the  static  true 
stress  -  true  strain  curve. 

A  correction  made  to  assumption  3  is  that  the  total  damaging  energy 
required  for  fatigue  fracture  is  not  constant  but  related  to  the  strain 
amplitude  per  cycle  as  follows: 


where  m  is  the  hysteresis  energy  exponent  which  Halford  and  Morrow  (Ref 
49)  derived  from  Griffith's  crack  theory.  Halford  and  Morrow's  value, 
m  =  4n,  is  accurate  for  4340  steel,  but  not  necessarily  for  other  metals. 
There  does  not  appear  to  be  a  relationship  between  strain  amplitude  and 
hysteresis  energy  to  failure  that  is  independent  of  material  properties. 

The  lack  of  a  relationship  between  strain  amplitude  and  hysteresis  energy 
does  not  invalidate  Equation  16  or  assumption  3,  it  simply  requires  experi¬ 
mental  determination  of  m  for  each  material  investigated.  Combining 
Equations  15  and  16  results  in  a  relationship  between  plastic  strain  and 
fatigue  life: 
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where  k  =  1/C  and  C,  the  fatigue-ductility  exponent  introduced  in  the 
Manson-Cof f in  law,  is  shown  to  be  l/(5n+l).  Again,  this  value  is  accurate 
only  for  4340  steel.  From  Radhakishnan  (Ref  28),  m  is  not  equal  to  4n, 
but  is  given  by  the  formula: 

m  =  k  -  (n  +  1) 

Both  n  and  C  (and  thus  k)  must  be  determined  experimentally.  Substituting 
£pu >  e  ,  1/2,  and  for  e^,  e^,  N^,  and  N^,  respectively,  and  rearranging 
Equation  15: 


N  =  1/2  (E  /E  ) 

f  p  pu 


produces  a  power  law  relation  between  fatigue  life  and  the  average  plastic 
strain  per  cycle.  In  Equation  17,  e^u  is  the  static  plastic  strain  fracture 
value  at  1/2  cycle.  Figure  7  demonstrates  the  relationship  between  fatigue 
life  and  average  plastic  strain  per  cycle  as  modeled  by  Equation  17. 

Criticism  has  been  levied  against  the  use  of  hysteresis  energy  as  a 
criterion  for  fatigue  failure.  The  first  criticism  is  that  attempts  to 
correlate  fatigue  life  with  tensile  properties  have  been  unsuccessful. 
However,  there  are  investigations  that  indicate  a  relationship  between 
the  fatigue  strength  and  parameters  of  the  true  stress-true  strain  curve 
(Ref  50  and  51).  The  second  source  of  criticism  is  due  to  the  difference 
between  monotonic  and  cyclic  stress -strain  curves.  This  is  a  minor  argu¬ 
ment,  since  the  difference  is  small  enough  to  neglect  when  the  two  curves 
are  compared  in  the  fatigue  stress  range  (Ref  28). 


Material  Properties.  Crack  growth  is  a  material-dependent  phenomenon. 
The  total  hysteresis  energy  to  failure,  W^,  and  plastic  energy  absorbed 
per  cycle,  dw,  were  presented  in  Equation  2.  For  most  materials,  the 
hysteresis  loop  is  relatively  stable  after  a  certain  initial  life.  This 
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stability  is  dependent  on  the  stress  level.  There  is  more  stability  at 
low  stress  levels  than  at  high  stress  levels.  In  some  cases,  like  carbon 
steels,  no  stable  loop  is  obtained  throughout  the  life  of  the  specimen, 
especially  when  the  stress  level  is  near  the  yield  strength.  A  consid¬ 
erable  amount  of  softening  is  observed  during  the  initial  stages,  followed 
by  a  hardening  process,  which  is  rapid  when  the  stress  level  is  high  and 
slow  when  the  stress  level  is  low.  Hence,  the  value  of  dw  per  cycle 
changes  with  stress  cycles  depending  on  the  hardening  and  softening  of 
the  material. 

The  empirical  relationships  and  mathematical  models  developed  to 
correlate  fatigue  life  with  mechanical  properties  of  metal  have  been 
criticized  for  depending  too  greatly  on  the  idealization  of  the  pro¬ 
perties  of  metals.  At  the  microscopic  level  where  fatigue  damage 
initiates,  all  metals  are  heterogeneous.  The  fatigue  resistance  of  a 
small  volume  of  material  will  differ  from  that  of  another  because  of 
inclusions,  differences  in  grain  size,  anisotropy,  orientation,  micro¬ 
residual  stress,  etc.  The  fatigue  behavior  of  a  real  material  is  the 
integrated  phenomenon  of  these  various  microscopic  factors.  The  dis¬ 
tribution  of  these  factors,  and  thus  the  metal’s  mechanical  properties, 
are  random  throughout  the  material.  Because  of  this  probabilistic 
distribution,  the  fatigue  limit  and  fatigue  life  at  a  given  stress  will 
vary  from  specimen  to  specimen. 

Plastic  Strain  and  Dislocation  Density 

An  empirical  relationship  between  plastic  strain  and  dislocation 
density  is  presented.  This  section  couples  together  fracture  criterion, 
based  on  plastic  work,  with  internal  damping,  based  on  dislocation 
density. 

Dislocation  Model  for  Stress-Strain  Behavior.  A  number  of  authors 
have  expressed  stress-strain  behavior  via  a  dislocation  model.  For  alpha 
iron  (o-Fe),  during  constant  strain-rate  tensile  experiments,  modeling 
the  yield-point  behavior  has  produced  a  variety  of  theories.  Cottrell 
(Ref  52)  explained  the  upper  yield  point  by  introducing  the  pinned  dis¬ 
location  concept  A  relatively  high  initial  stress  was  required  to  unpin 
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dislocations,  but  once  unpinned  they  could  move  at  a  lower  stress  level 
In  his  approach,  Hahn  (Ref  53)  explained  the  yield-point  phenomenon  as  a 
result  of  the  rapid  increase  of  mobile  dislocations  when  yielding  begins 
A  model  developed  somewhat  later  by  Bergstrom  (Ref  54),  emphasizes  the 
A-e  relationship  and  is  based  on  dislocation  behavior  during  deformation 

Bergstrom  has  five  assumptions  for  his  model: 

1.  The  true  flow  stress-dislocation  density  relation  is  given  by: 

o  =  o.+  ctGbA1/2  (18) 

10 


where  a  and  G  are  constants.  A  review  of  experimental  data  in  support 
of  Equation  17  is  provided  by  Otte  and  Hren  (Ref  55). 

2.  Both  mobile  and  immobile  dislocations  are  present  at  any  true 
strain,  e. 

3.  The  mobile  density,  L,  is  strain  independent  and  much  smaller 
than  the  immobile  density,  A^. 

4.  The  variation  of  A  with  €  is  determined  by  the  creation,  the 
immobi 1 izat ion,  the  remobilization,  and  the  annihilation  of  dislocations 

5.  Those  effects  that  arise  from  changes  to  the  scale  or  arrange¬ 
ment  of  the  dislocation  structure  and  that  do  not  affect  assumption  4 
above,  are  negligible. 


Bergstrom  derived  the  following  differential  equation,  giving  the 
dislocation  density,  A,  as  a  function  of  plastic  strain,  e  : 


dA 

de 


U  -  XA 


(19) 


where  U  is  the  rate  of  immobilization  of  mobile  dislocations  and  x  is 
the  stra in- independent  probability  for  the  remobilization  of  immobile 
dislocations,  and: 
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1/ibs 


U.  = 
1 


where  #  is  an  orientation  factor  (0.5  and  0.32  for  h.c.c.  and  f.c.c., 
respectively)  and  s  is  the  dislocation  mean  free  -path.  Bergstrom  (Ref 
54)  contends  that  a  strain- independent  s  is  a  good  approximation  for 
b.c.c.  metals  if  X  is  small  enough  for  the  product  xA  to  be  negligible 
at  small  strains.  With  a  strain-independent  U.,  the  integration  of 


Equation  19  results  in: 


A  =  [1  -  exp( -xe  ) ]  +  AQexp(-XEp) 


with  A  =  A.  at  t  =0. 

0  p 

As  a  modification  to  the  Bergstrom  theory,  Vetter  and  van  den  Beukel 
(Ref  56)  assume  the  dislocation  mean  free  patli  is  of  the  order  of  and 
proportional  to  the  dislocation  cell  size,  d.  The  cell  size  decreases 
with  increasing  strain,  and  has  been  shown  to  depend  on  the  dislocation 
density  according  to: 

d  =  ± 


where  K  =  20  for  iron.  Therefore: 


s  =  f 


where  f  is  a  constant  of  order  1.  Inserting  Equation  20  into  Equation 
19  and  integrating 


A1/2  =  ~  [1  -  exp(-l/2xtp)  1  +  A0^2exp(-l/2XEp) 


where  A'  =  l/$bKf. 
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Damping  Mechanism-Fatigue  Criterion.  Equation  21  provides  a 

relationship  between  plastic  strain  amplitude  and  dislocation  density 

and  is  the  link  to  combining  the  damping  mechanism  (Equation  13)  with 

the  fatigue  criterion  (Equation  17).  Rearranging  Equation  21  so  as  to 

erpress  e  as  a  function  of  A, 

P 


-t1/2  X  -  A* 

V'2*  -  *• 


For  the  fatigue  criterion,  as  described  by  Equation  17,  the  average 
plastic  strain,  e  ,  is  the  sum  of  the  plastic  strain  per  cycle,  e  , 
divided  by  the  number  of  elapsed  stress  cycles.  In  equation  form: 


where  N  is  the  number  of  elapsed  cycles.  The  right-hand  side  of  Equation 
22  is  then  substituted  for  e  in  Equation  23  and  the  fatigue  criterion 
becomes  a  function  of  dislocation  density: 


l  (r 


/i  ,  .  r 

i  X  -  A 
A0'/2  X  -  A' 


Damping  was  previously  shown  to  be  dependent  on  dislocation  density 

(see  Equation  13);  however,  this  dependency  consists  of  two  dislocation 

density  terms:  A  and  A  .  It  is  now  assumed  that  the  plastic  zone  dis- 
e  p 

location  density,  A^,  is  responsible  for  all  plastic  strain  directly 
related  to  crack  growth,  so  that  no  crack  growth  exists  in  the  elastic 
region.  That  is  to  say.  Equation  22  may  be  rewritten  as: 


,A 

In  M 


Ac/'2  *  '  A' 
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Equation  25  is  consistent  with  the  fatigue  process  N  stage  (Ref  11,  12) 
where  voids  form  in  the  plastic  zone  and  coalesce  into  microcracks. 

This  situation  is  illustrated  in  Figure  8  and  is  modeled  by  Ishikawa 
(Ref  57)  as  a  slow  crack  extension  within  the  plastic  zone,  where  the 
plastic  zone  is  a  uniform  distribution  of  dislocations.  Ishikawa  derives 
the  plastic  strain  in  the  plastic  zone  from  the  crack  extension  force 
as : 


e  =  «(s/r  )  (26) 

P  P 

where  a  is  approximate  unity. 

Rearranging  Equation  13  into: 

AC  V  f  E2dV  -  A  [V  f  e28^(e )dV  +  C^V  f  £2dVl 

- - 1 - _e - b - -  (27) 

V  fp  E20jj(£)dV  +  C3Vp  f  £2dV 

solves  for  the  plastic  zone  dislocation  density  A  in  terms  of  the  cor- 

— c  ^ 

rected  total  decrement,  A  ;  the  elastic  region  dislocation  density,  A^; 
the  plastic  zone  radius,  r  ;  and  the  strain  amplitude,  e.  Dislocation 
density  increases  for  some  metals  at  stresses  below  the  macroscopic  yield 
(Ref  58)  and  under  these  conditions  it  is  not  possible  to  eliminate  A 

e 

as  a  variable.  However,  for  metals  whose  dislocation  density  does  not 
increase  appreciably  until  yield,  A^  may  be  used  in  place  of  A^  (Ref 
56). 

Equation  17  presents  fatigue  life  as  a  function  of  average  plastic 
strain.  Equation  24  determines  average  plastic  strain  from  dislocation 
density  in  the  crack  tip  plastic  zone.  Equation  27  gives  plastic  zone 
dislocation  density  from  internal  damping.  With  Equations  17,  24,  and 
27,  the  objective  of  developing  a  theoretical  relation  between  fatigue 
and  damping  has  been  accomplished.  The  next  requirement  is  to  test  the 
validity  of  Equations  17,  24  and  27. 


A 

P 
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EXPERIMENTAL  APPROACH 


Experimental  verification  of  the  theoretical  development  of  a 
relationship  between  crack  growth  and  internal  damping  begins  with  the 
determination  of  material  properties,  including  damping,  and  concludes 
with  fatigue  testing  of  1018  and  4340  steel  beam  specimens. 

Purpose 

The  main  purpose  of  experiments  is  to  test  hypotheses.  The  fatigue 
tests  used  in  this  research  will  determine  whether  internal  damping 
changes  with  crack  growth  as  suggested  by  Equations  17,  24,  and  27. 

Damp i ng  Mechanism-Fatigue  Criterion  Relationship.  Successful  ver¬ 
ification  of  Equations  17,  24,  and  27  in  a  test  series  of  one  material 
and  specimen  type  is  no  assurance  that  the  proposed  constitutive  rela¬ 
tionship  between  crack  propagation  and  internal  damping  is  valid  for  all 
material  and  specimen  types.  It  merely  states  that  the  relationship 
between  fatigue  and  damping  is  valid  for  the  material  and  strain  dis¬ 
tribution  tested.  A  measure  of  validity,  rather  than  complete  verifi¬ 
cation,  is  required  at  this  point.  Determining  whether  or  not  the 
relationship  can  be  established  for  a  particular  material  and  specimen 
type  is  the  first  step.  This  relationship  was  built  through  a  process 
of  steps  culminating  in  Equations  24,  25,  and  27.  As  many  of  these  steps 
as  possible  must  be  independently  verified  to  assist  in  determining  the 
problem  areas  if  discrepancies  arise.  Therefore,  the  primary  purpose  of 
testing  is  to  determine: 

1.  Ts  plastic  work  related  to  fracture  as  suggested  by  Equation 
15? 

2.  Is  plastic  strain  related  to  fracture  by  Equation  17? 

3.  Does  Equation  13  accurately  model  internal  damping  as  a 
function  of  strain  and  plastic  zone  and  elastic  region 
dislocation  densities? 
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4.  Will  the  Bergstrom  model  relate  plastic  strain  and  dislocation 
density  as  expressed  in  Equation  22? 

5.  Do  all  the  previous  steps  combine  to  relate  crack  propagation 
to  internal  damping? 

The  five  questions  listed  above  interrelate  internal  damping,  crack 
propagation,  plastic  strain,  and  dislocation  density.  These  four  param¬ 
eters  must  all  be  evaluated  during  testing  in  order  to  verify  the  internal 
damping- fatigue  criterion  relationship. 


Fatigue  Life,  Damage  Assessment,  and  Fracture  Prediction.  A  step 
beyond  verifying  the  internal  damping-crack  growth  relationship  is  to 
predict  fatigue  specimen  life  and  the  subsequent  specimen  failure.  The 
fracture  prediction  process  entails  extrapolating  the  average  plastic 
strain  calculated  from  the  damping  data  to  arrive  at  the  predicted  total 
number  of  fatigue  cycles  to  fracture,  ,  as  presented  in  Equation  17. 

A  simple  approach  to  damage  assessment  is  to  use  the  ratio  of  fatigue 
cycles  completed,  Nc,  to  fatigue  cycles  to  fracture,  N^,  as  predicted 
by  the  fatigue  life  estimate  of  Equation  17.  In  equation  form: 


D%  =  100%  (N  /N  ) 
c  fp 


(28) 


where  D%  is  the  damage  assessment  in  percent.  A  comparison  of  this 
estimate  with  the  actual  crack  growth  will  suffice  as  verification. 

This  linear  approach  to  damage  assessment  should  prove  satisfactory  for 
low  cycle  fatigue  using  constant  crack  growth  rates. 

The  linear  prediction  of  fatigue  damage  represented  by  Equation  28 
is  adequate  for  loads  of  a  constant  stress  intensity  factor;  however,  it 
is  inaccurate  for  other  types  of  loads.  The  life  of  a  structure  under¬ 
going  fatigue,  when  based  on  plastic  work,  is  given  by  Equation  17. 
Specifically,  for  a  varying  plastic  strain  per  fatigue  cycle,  t 
the  fraction  of  fatigue  life  expended  for  a  structure  is  given  by: 


pr 


LF  = 


(29) 
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where  e^u  is  the  static  plastic  strain  fracture  value  at  1/2  cycle.  As 
the  crack  growth  to  failure,  a^,  divided  by  the  crack  growth  in  one  cycle, 
Aa ^ ,  results  in  fatigue  cycles  to  failure,  N^,  substitution  in  Equation 
17  produces  the  following  relation: 


t 

PJ 


,  Aa , 

t  i  -J- 
pu  \af 


(30) 


Substitution  of  e  .  into  Equation  21  determines  the  plastic  zone  dis- 
PJ  th 

location  density  for  the  j  fatigue  cycle.  Combined  with  an  estimate 
of  the  plastic  zone  radius  from  Equation  11  or  12,  an  approximation  of 
the  internal  damping  can  be  made  from  Equation  13. 

Equation  29  is  presented  as  a  more  accurate  and  reliable  estimation 
of  fatigue  damage  than  Equation  28.  Successful  verification  of  Equation 
29  will  meet  the  objective  of  establishing  nondestructive  evaluation 
techniques  based  on  internal  damping. 

Test  Arrangements 

The  two  test  arrangements  described  here  are  designed  to  answer  the 
questions  raised  above.  In  so  doing,  these  experiments  will  verify  the 
proposed  relationship  between  fatigue  and  damping  for  1018  and  4340  steel. 

Fatigue  Test.  To  accomplish  the  desired  goals  of  relating  damping 
to  crack  growth  and  predicting  fatigue  life,  a  series  of  fatigue  tests 
is  required.  These  tests  monitor  crack  growth,  damping,  load,  strain, 
and  number  of  fatigue  cycles. 

The  fatigue  test  schematic  is  shown  in  Figure  9,  the  fatigue  beam 
specimen  in  Figure  10,  and  the  test  procedure  is  described  in  detail  in 
Appendix  A.  The  following  description  briefly  summarizes  the  fatigue 
test  procedure. 

An  MTS  test  machine  applies  a  cyclic  load  to  a  4340  beam  specimen 
(100  pounds  for  damping  datum  and  up  to  1,400  pounds  for  fatigue).  The 
beam  specimen  is  initially  notched  and  uncracked.  The  fatigue  load  has 
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a  stress  ratio  of  zero  and  is  continuously  adjusted  to  maintain  a  constant 
stress  intensity  factor  during  the  fatigue  process.  Crack  length  data 
is  collected  during  data  collection  cycles  from  crack  propagation  gages 
on  either  side  of  the  specimen.  The  crack  length  measurement  is  used  to 
update  the  fatigue  load  setting.  The  stress  intensity  factor  is  calculated 
from  (Ref  59): 


Kj  =  o/ia  F(a/b)  (31) 

2 

where  o  =  6M/b 

M  =  moment 

b  =  specimen  thickness 

and 

F(a/b)  =  1.122  -  1. 4(a/b)  +  7.33(a/b)2  -  13.08(a/b)3 

+  14.0(a/b)4  (32) 

is  an  empirical  formula  accurate  to  0.2  percent  for  a/b  <0.6.  In  the 

case  of  a  four-point  beam,  the  moment  is  given  by  M  =  1/2P1  ,  where  P  is 

d 

the  total  applied  load  and  1  is  the  distance  from  a  knife  edge  to  the 
nearest  load  point.  Substituting  for  0  in  Equation  31  and  solving  for 
P: 


(33) 


P  of  Equation  33  is  calculated  after  every  data  cycle  and  a  new  fatigue 
load  setting  is  communicated  to  the  MTS  test  machine. 

A  vibration  generator  attached  to  the  specimen  excites  the  beam  in 
a  range  of  5  to  20  Hz  around  the  fundamental  frequency,  w  j.  For  example, 
if  the  specimen  resonates  at  320  Hz,  a  frequency  range  from  310  to  330 
Hz  is  swept  by  the  vibration  generator.  While  the  beam  is  being  excited, 
strain  gage  data  are  input  to  a  spectrum  analyzer  that  converts  the  data 
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from  the  time  domain  to  the  frequency  domain  and  then  transfers  it  to 
the  Lab-Datax  computer  for  damping  calculations.  The  same  data  also 
provide  the  fundamental  frequency  of  the  beam.  For  each  damping  data 
cycle  the  following  data  are  collected: 

1.  Fatigue  cycle  number,  N 

2.  Fatigue  load,  P 

3.  Maximum  strain,  e 

a 

4.  Minimum  strain,  e 

min 

5.  Crack  length,  a 

6.  Stress  intensity  factor, 

7.  Damping  decrement,  A 

8.  Transfer  function  amplitude  ratio  at  u  ,,  TF 

nl 

The  data  are  saved  in  a  file  on  disk  to  be  transferred  at  a  later  time 
to  a  mainframe  for  processing  and  plotting. 

All  fatigue  tests  will  begin  with  a  2000  load  cycle  damping  datum 
phase.  During  this  phase,  approximately  50  damping  data  cycles  will  be 
performed  and  the  resultant  average  damping  value  will  be  interpreted  as 
the  normal  damping  state  of  the  uncracked  specimen  as  configured  and 
mounted  in  the  fixture.  The  load  during  this  phase  will  be  approximately 
100  pounds.  The  light  loading  is  sufficient  to  detect  minor  slippage  at 
the  knife  edges  and  make  appropriate  adjustments  to  the  clamping  bolts. 

Fatigue  tests  of  1018  steel  beam  specimens  are  similar  to  those 
described  for  the  4340  steel  specimens.  The  major  difference  is  in  the 
fatigue  load.  For  the  1018  steel  specimens,  the  fatigue  load  is  approx¬ 
imately  800  pounds.  Also,  the  load  is  held  constant.  Therefore,  the 
stress  intensity  factor  increases  as  the  crack  grows.  The  stress  inten¬ 
sity  factor  is  held  constant  for  the  4340  steel  specimens. 

Microcracks  initiate  in  the  vicinity  of  a  stress  concentration  as  a 
consequence  of  high  load  stress.  High  load  stress  forces  dislocations 
to  pile  up  on  a  slip  plane.  A  microcrack  is  nucleated  as  a  result  of 
the  short-range  interactions  of  a  few  dislocations  at  the  head  of  the 
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pile.  The  remaining  dislocations  are  presumed  to  produce  the  high  stress 
concentration  needed  to  force  the  leading  dislocations  to  form  a  micro¬ 
crack  (Ref  11,  60,  61,  62).  In  ductile  metals,  where  dislocations  are 
easily  generated,  crack  nuclei  continue  to  form  with  increasing  amounts 
of  strain,  and  these  nuclei  grow  slowly  into  microcracks  that  eventually 
join  to  form  a  macroscopic  crack.  Elasto-plastic  analysis  of  cracks 
often  ignore  the  microcrack  phase  and  directly  model  crack  propagation 
with  dislocation  strips  or  arrays  (Ref  57,  63,  64). 

For  a  cyclic  load  having  a  compressive  mean  stress,  microcracks  are 
prevented  from  opening,  retarding  development,  and  impeding  the  coalescence 
of  microcracks  into  macrocracks.  On  the  other  hand,  a  tensile  mean  stress 
tends  to  open  a  microcrack  and  spur  crack  development  (Ref  60).  The 
fatigue  tests  of  4340  steel  beams  are  designed  to  generate  a  tensile 
mean  stress.  The  microcracks  developed  in  these  specimens  remain  open, 
under  a  zero  stress  condition,  due  to  residual  plastic  stress. 

Within  the  plastic  region,  elliptical  microcracks  have  displacement 
in  two  modes:  opening  and  shear.  Opening  displacement  spreads  and  then 
closes  the  microcrack.  This  mode  creates  internal  damping  through  energy 
loss  when  the  surfaces  are  in  contact.  Shear  displacement  moves  one 
microcrack  surface  parallel  to  the  other.  Damping  is  incurred  through 
the  dissipation  of  heat  as  the  two  surfaces  rub  past  one  another. 

The  following  assumptions  were  set  forth  by  Budiansky  and  O'Connell 
(Ref  65),  in  their  analysis  of  randomly  distributed  cracks: 

1.  Microcrack  closure  effects  are  ignored. 

2.  Microcracks  have  small  openings  between  their  opposite  faces. 

3.  Microcrack  edges  are  considered  blunt. 

4.  Small  stresses  do  not  produce  contact  between  microcrack  faces. 

5.  Macroscopic  incremental  stress-strain  relationship  is  linear. 

6.  General  elliptic  pi  an  form. 
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Having  a  mean  tension  fatigue  stress,  microcracks  in  the  4340  beam 
specimens  meet  the  condition  of  assumption  2.  Given  the  small  damping 
data  cyclic  stress  (on  the  order  of  500  psi)  applied  to  the  specimen 
(described  in  Appendix  A),  it  is  reasonable  to  presume  that  condition  4 
is  also  met.  Meeting  conditions  2  and  4  suggest  that  the  microcrack 
surfaces  do  not  touch  and  thereby  dissipate  energy.  Without  dissipation 
of  energy,  the  microcracks  cannot  contribute  to  internal  damping. 

Two  alterations  to  the  fatigue  test  arrangement  of  Appendix  A  could 
result  in  microcracks  contributing  to  internal  damping.  First,  the  mean 
fatigue  stress  would  have  to  be  in  a  range  that  would  allow  the  microcrack 
surfaces  to  move  relative  to  one  another  under  a  small  data  collection 
cyclic  loading.  The  motion  could  be  shear  mode,  opening  mode,  or  a  com¬ 
bination  of  the  two.  Second,  the  data  collection  cyclic  loading  would 
have  to  be  large  enough  to  overcome  the  residual  plastic  strain  holding 
the  microcracks  open  and  cause  the  microcrack  surfaces  to  rub  against 
one  another. 

Internal  Damping  Test.  Material  parameters,  required  to  calculate 
dislocation  damping  in  the  plastic  zone  from  Equation  27,  are  determined 
experimentally  by  internal  damping  tests.  This  section  summarizes  the 
test  arrangement  for  determining  internal  damping  material  parameters. 

Physical  properties  of  materials  are  determined  by  repeatable  exper¬ 
imental  testing.  Internal  damping,  being  a  physical  property,  must  be 
determined  by  a  series  of  tests  for  each  material  of  interest.  To  answer 
the  question  whether  Equation  13  accurately  models  internal  damping  as  a 
function  of  strain  and  plastic  zone  and  elastic  region  dislocation  densi¬ 
ties,  testing  of  4340  steel  at  various  levels  of  strain  is  required. 

Figure  ll  is  a  schematic  of  the  damping  test  arrangement  and  Figure 
12  shows  the  damping  specimen.  Appendix  B  contains  a  description  of  the 
test  procedure.  The  following  description  summarizes  this  procedure.  A 
sine  load  function  is  applied  to  the  symmetric  cantilever  beam  specimen 
furnishing  a  stress  ratio  of  negative  one.  An  accelerometer  measures 
the  load  function  while  strain  gages  on  the  specimen  measure  the  beam 
response.  The  ratio  of  strain  to  acceleration  at  the  resonant  frequency 
of  the  beam  provides  a  transfer  function  ratio.  From  this  ratio,  internal 
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damping  for  the  material  may  be  plotted  against  strain  levels,  thus 
obtaining  plots  similar  to  Figures  2  and  4.  The  data  collected  from 
these  tests  will  include  decrement,  A,  and  strain,  e.  From  data  reduc¬ 
tion,  the  parameters  C^,  C^,  and  of  Equations  14  and  26  are  estimated 
for  later  application  in  refining  the  fatigue  test  data  reduction. 

Material  Properties 

Frior  to  performing  fatigue  and  damping  tests,  testing  was  done  to 
determine  the  mechanical  properties  of  4340  steel.  The  chemical  composi¬ 
tion,  heat  treatment,  and  engineering  tensile  properties  of  the  SAE  4340 
steel  used  in  this  investigation  are  presented  in  Table  1.  The  SAE  4340 
steel  was  selected  as  the  test  material  for  several  reasons:  (1)  a  sig¬ 
nificant  amount  of  crack  growth  data  had  been  previously  collected,  (2) 
a  small  plastic  zone  on  the  order  of  0.01  to  0.03  inch  developed  at  the 
notch,  and  (3)  there  was  continued  usage  of  the  material  in  applications 
with  conditions  leading  to  cyclic  fatigue.  The  large  amount  of  crack 
growth  data  available  for  4340  steel  allows  quick  verification  of  fatigue 
test  results.  Conversely,  1018  steel  produces  a  large  plastic  zone  that 
sharply  contrasts  with  4340  steel. 

Testing  and  subsequent  calculations  (Ref  66  and  48)  were  done  to 
determine  the  true  stress-true  strain  relationship  of  the  materials. 
Figure  13  depicts  this  relationship  and  has  a  strong  resemblance  to 
Figure  6,  the  true  stress-strain  curve  of  4340  steel  used  by  Feltner  and 
Morrow  (Ref  47).  The  true  stress-plastic  strain  curve  representative  of 
Figure  13  may  be  expressed  as: 


0 
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k  e 
P  P 


(34) 


where  the  strength  coefficient  k^  =  230.3x10'  ps i ,  and  the  strain  hardening 
exponent  n  ~  0.0902.  Coefficients  k^  and  n  were  computed  using  a  least 
squares  technique  (Ref  67). 

Inserting  Equation  21  into  Equation  18  produces  a  second  equation 

relating  stress,  o  and  plastic  strain,  e  , 

P 
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=  0(0)  +  —  aGb[ 1  -  exp(-l/2XE  )] 

x  p 


where  o(0)  is  the  stress  at  zero  plastic  strain.  A  simplex  algorithm 

(Ref  68)  was  used  to  fit  the  data  to  a  curve  which  resulted  in  A'  = 

6  "  1  3 

6.128x10  in.  ,  o(0)  =  162.8x10  psi,  and  X  =  8.613  assuming  a  =  0.8 

_  g 

and  the  Burger's  vector,  b  =  10  .  The  shear  modulus  is  calculated  from 

(Ref  69): 


G  = 


E 

2( 1+v) 


(36) 


where  the  elastic  modulus  has  been  experimentally  determined  as 
E  =  28.0x10^  psi  and  Poisson's  ratio  is  v  =  0.29  (Ref  66).  Equation  36 
produces  G  =  10.9x10^  psi.  From  Feltner  and  Morrow  (Ref  47)  the  total 
hysteresis  energy  for  static  fracture,  based  on  Equation  34,  is  given 
by: 
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,  n+1 

k  e 
P  Pu 
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(37) 


or  based  on  Equation  35, 

A '  2A ' 

W  =  [o(0)  +  aGb]  e  +  aGb[exp(-l/2XE  )  -  1]  (38) 

u  X  pu  ^2  1  pu  ' 

The  static  test  results  are  summarized  in  Table  2. 

The  final  preliminary  tests  determined  the  critical  stress 

intensity  factor  of  4340  steel,  K  .  Three  compact  tension  specimens 

*c  3  1/2 

(Ref  70)  were  tested,  producing  a  =  76.5x10  psi-in.  .  From 
Equation  12  and  the  yield  stress  given  in  Table  1,  the  plastic  zone  r 
would  be  expected  to  reach  a  maximum  of  0.016  inches  for  a  specimen 
experiencing  plane  strain  conditions.  Equation  27  integrates  the  dis¬ 
location  density  and  strain  functions  over  volumes  of  the  plastic  and 
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elastic  beam  material.  With  an  elastic  volume  of  about  5.6  in.  versus 

*6  3 

a  plastic  volume  of  about  800  x  10  in.  ,  the  ratio  of  plastic  to  elastic 
material  is  on  the  order  of  150  x  10  This  ratio  suggests  a  large 
plastic  zone  dislocation  density  or  strain  level  is  required  to  produce 
a  measurable  internal  damping  change. 

These  values,  together  with  the  results  from  the  internal  damping 
tests,  will  be  used  in  the  mathematical  model  to  verify  Equations  17, 

24,  and  27. 


TEST  RESULTS 


The  results  of  the  four-point  beam  fatigue  tests,  and  brief  results 
of  the  internal  damping  tests,  are  presented  here.  Detailed  results  of 
the  internal  damping  tests  are  in  Appendix  E. 

In  addition  to  the  dislocation  damping  discussed  previously,  an 
additional  damping  mechanism  is  analyzed  for  the  4340  steel  beam  speci¬ 
mens.  This  mechanism  dominates  the  damping  response  observed  in  these 
specimens . 

Internal  Damping  Tests 


To  validate  the  fatigue  test  results  and  mathematical  model, 
internal  damping  tests  were  conducted.  This  section  summarizes  the 
results  described  in  detail  in  Appendix  E,  and  is  based  on  the  test 
procedure  presented  in  Appendix  B. 

The  coefficients  associated  with  dislocation  damping,  and  C^, 

are  used  in  Equations  13  and  14,  and  subsequently  in  Equation  27. 
Modeling  of  these  equations  is  established  in  Appendix  C,  step  3.8. 
Appendix  C  is  an  algorithm  for  simulating  internal  damping  in  prepara¬ 


tion  for  the  mathematical  model. 


Estimates  of  material  constants  are  C  =  54  x  10  and  C  =  144  x 
*6  i  z 

10  Based  on  restrictions  of  the  Granato-Liicke  equation  (Equation  5), 


the  maximum  beam  strain  during  damping  data  collection  is  limited  to  36 
x  10  that  is,  one  fourth  of  C^.  From  Figure  14,  for  a  simple-supported 


beam  vibrating  in  the  fundamental  mode,  the  damping  ratio  for  the  maximum 
strain  is  172  x  10  ^ .  This  small  amount  of  damping  is  difficult  to  mea¬ 
sure  accurately,  a  point  that  becomes  more  apparent  in  the  next  section, 
which  documents  the  analysis  of  the  four-point  beam  fatigue  tests. 

Fatigue  Tests 

The  two  materials  used  in  the  fatigue  tests  were  1018  and  4340 
steel.  The  1018  steel  beam  specimens  were  fatigued  under  a  constant 
load  and  an  increasing  stress  intensity  factor.  The  4340  steel  beam 
specimens  were  fatigued  under  a  decreasing  load  and  a  constant  stress 
intensity  factor. 

4340  Steel  Beam  Fatigue  Plots.  Four  beam  specimens  were  fatigue 
tested  in  a  four-point  load  arrangement.  Figures  15  through  18  show  the 
plotted  crack  extension,  resonant  frequency,  and  damping  ratio  data  for 
each  of  the  specimens.  The  figures  share  common  peaks  in  damping  and 
fluctuations  of  resonant  frequency  at  approximately  300,  270,  and  245 
Hz.  Since  the  tests  were  automated,  no  observations  are  available  to 
explain  these  phenomena  in  the  data.  Due  to  differences  in  the  initial 
tightness  of  the  clamping  bolts  at  the  knife  edges,  the  beams  began  their 
tests  at  different  resonant  frequencies.  The  specimen  crack  lengths 
also  varied  from  beam  to  beam  at  these  frequencies.  Because  of  this,  it 
does  not  appear  that  the  fluctuations  were  due  to  the  specimens  them¬ 
selves.  One  possible  explanation  is  that  some  part  of  the  test  equip¬ 
ment  or  fixture  on  the  MTS  test  machine  resonates  at  these  frequencies. 

All  four  plots  indicate  a  definite  increase  in  damping  as  the  crack 
grows.  This  upward  trend  in  damping  is  consistent  with  the  increase  in 
energy  losses  in  metals  with  propagating  cracks  as  measured  by  Charles, 
Appl,  and  Francis  (Ref  33).  Freudenthal  (Ref  71),  noting  the  sharp 
increase  in  damping  and  fracture  that  ensued,  offered  the  explanation 
that  the  increase  is  an  expression  of  the  formation  and  spreading  of 
macroscopic  cracks  under  considerable  local  deformation. 
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Specimen  A.  Crack  extension,  resonant  frequency,  and  damping 
ratio  data  are  presented  in  Figure  15  for  this  first  specimen  tested. 
Fatigue  loading  was  set  at  1250  pounds  with  a  resulting  K  g  =  26.1  and 
SD  =  2. 16  ksi-in. where  SD  is  the  standard  deviation.  Due  to  a  pro¬ 
gramming  error  in  the  calculation  of  Equation  32,  the  stress  intensity 
factor  was  not  held  constant.  A  decision  was  made  to  fail  the  specimen 
while  the  crack  growth  was  still  being  tracked  by  the  crack  propagation 
gages  in  order  to  assess  their  accuracy.  Table  3  summarizes  the  data  at 
failure.  The  crack  propagation  gage  recorded  a  crack  depth  of  approxi¬ 
mately  0.25  inch  at  failure,  indicating  an  error  on  the  order  of  10  per¬ 
cent.  This  error  is  probably  associated  with  the  error  in  measuring  the 
location  of  the  crack  gages  on  the  specimen  prior  to  the  start  of  the 
test . 


Specimen  B.  The  damping  datum  was  established  at  a  damping 

ratio  of  460  x  10  ^  and  the  load  was  initialized  at  1250  pounds.  Fatigue 

test  data  are  shown  in  Figure  16.  The  average  stress  intensity  factor 

1/2 

was  34.4  with  SD  =  1.37  ksi-in.  .  The  average  crack  growth  rate, 
(da/dN),  over  the  width  of  the  crack  propagation  gage  was  14.9  x  10 

Specimen  C.  Figure  17  contains  the  crack  extension,  resonant 

frequency,  and  damping  ratio  data.  The  damping  datum  was  established  at 

a  damping  ratio  of  460  x  10  ^  over  50  data  cycles.  The  fatigue  load  was 

1/2 

initialized  at  1250  pounds.  K  =  30.5  and  SD  =  1.91  ksi-in.  with 

ave 

the  crack  growth  rate,  (da/dN),  at  11.9  x  10 


Specimen  D.  The  data  from  the  last  test  are  presented  in 

Figure  18.  The  damping  was  established  at  a  damping  ratio  of  460  x  10  ^ 

and  the  load  was  initialized  at  1400  pounds.  An  average  stress  intensity 

1/2 

factor  of  36.3  with  a  standard  deviation  of  1.48  ksi-in.  was  produced. 
Over  the  width  of  the  crack  gage  the  crack  propagation  rate  was  17.0  x 
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Past  experience  with  crack  propagation  gages  had  produced  inconsis¬ 
tent  results.  The  adhesive  selected  to  glue  the  gages  to  a  specimen 
often  failed,  resulting  in  erratic  crack  length  data.  The  adhesive  used 
on  the  4340  steel  beam  specimens  proved  very  reliable  and  secured  crack 
length  data  within  100^  of  an  inch. 

Analysis  of  Damping  Results  in  4340  Fatigue  Specimens.  Results 
from  the  damping  tests  are  applied  to  the  fatigue  tests  through  the  esti¬ 
mates  of  the  material-dependent  constants  of  Equation  E-3  in  Appendix  E 
and  usage  of  the  frequency-dependent  damping  equation. 

The  plastic  strain  in  the  plastic  zone  at  the  fatigue  crack  fip  was 
estimated  at  400  x  10  ^  in. /in.  for  specimens  B  and  C  and  500  x  10  °  for 
specimen  D.  This  estimate  is  based  on  Equation  17  and  assumes  a  constant 
plastic  zone  size.  For  fatigue  specimens  B,  C,  and  D  the  assumption  of 
a  constant  plastic  zone  size  is  valid  for  the  period  of  testing  when  the 
crack  length  was  within  the  crack  propagation  gage  range.  Once  outside 
of  this  range,  erroneous  crack  length  data  were  received  by  the  control 
program  and  the  applied  fatigue  load  did  not  preserve  a  constant  stress 
intensity  factor.  This  small  plastic  strain,  together  with  the  small 
plastic  zone  radius  (on  the  order  of  0.003  inch)  suggests  a  small  con¬ 
tribution  (on  the  order  of  10  to  the  overall  internal  damping  by  the 
dislocations  in  the  plastic  zone. 

Figures  15  through  18  clearly  indicate  a  significant  increase  in 
the  internal  damping  as  each  test  progressed.  However,  the  cause  of 
this  increase  cannot  be  associated  with  the  dislocation  density  of  the 
plastic  zone  at  the  fatigue  crack  as  hypothesized.  Another  mechanism 
must  be  found  to  explain  this  phenomenon. 

Alternate  Damping  Mechanism.  The  damping  behavior  observed  in 
Figures  15  through  18  is  not  consistent  with  phenomenon  associated  with 
dislocation  density  in  a  plastic  zone.  For  a  constant  plastic  zone  radius, 
this  damping  should  be  relatively  constant.  Damping  in  all  specimens 
increased  with  increasing  crack  length,  suggesting  that  damping  is  related 
to  the  crack  length  itself.  Plots  of  damping  versus  maximum  beam  slope 
at  the  crack  (Figures  19  through  21)  reveal  a  straight  line.  Although 
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the  slope  is  different  for  the  three  figures,  it  does  suggest  that  the 
cause  of  the  damping  is  related  to  the  crack  length  or  rotation  of  the 
cracked  surface.  A  relation  based  on  rotation  of  the  beam  at  the  crack 
surface  has  the  form: 


<t>'  (1/2) 

max 


An  initial  conjecture  is  that  the  crack  surfaces  rub  together  when  the 
beam  vibrates.  Figure  22  models  the  damping  versus  beam  slope  relation. 
This  model  suggests  the  following  boundary  equation  using  a  constant, 
c  ,  ,  multiplying  the  beam  slope  velocity  to  represent  the  crack  surface 
dissipation  function: 


F.Iw"(  1/2)  =  cck(a)  w'(l/2)  +  kja)  w’(l/2)  (39) 

This  is  an  empirical  relation,  adequate  for  the  model  in  this  report, 
but  requiring  additional  analysis  at  some  later  date. 

1018  Steel  Beam  Fatigue  Test  Results.  Figure  23  is  a  plot  of  data 
collected  from  a  1018  steel  beam  subjected  to  fatigue  loading.  The 
beam  specimen  has  a  shallow  increase  in  damping  between  17,500  to  20,000 
cycles;  however,  it  is  not  clear  if  this  change  in  damping  can  be  attrib¬ 
uted  to  the  fatigue  phenomenon  or  if  it  is  caused  by  an  external  effect. 
One  possibility  is  that  this  additional  damping  is  due  to  the  crack 
motion,  as  reflected  in  Equation  39  and  the  coefficient  c  A  scenario 
for  this  would  be  that  the  plastic  zone  is  initially  small  and  the  crack 
surfaces  rub  as  the  beam  is  vibrated  to  collect  damping  data.  As  the 
crack  grows  and  the  stress  intensity  factor  increases,  the  beam  forms  a 
plastic  hinge  at  its  midpoint.  The  plastic  hinge  permanently  opens  the 
crack,  not  allowing  the  surfaces  to  rub  and  effectively  lowers  the  damp¬ 
ing  ratio.  Excluding  the  damping  increase  mentioned  above,  the  overall 
damping  ratio  remains  roughly  constant  up  to  25,000  cycles,  whereupon  it 
increases  sharply. 
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Test  results  of  1018  and  4340  steel  specimens  demonstrate  two  dif¬ 
ferent  damping  phenomenon.  Results  from  the  1018  steel  fatigue  tests 
are  consistent  with  the  analytical  approach  presented  previously  in  this 
report.  That  is,  the  damping  is  relatively  stable  for  a  small  plastic 
zone,  and  increases  rapidly  with  an  expanding  plastic  zone.  To  validate 
the  1018  steel  results  and  analyze  the  4340  steel  results,  the  specimens 
must  be  modeled  and  crack  growth  and  damping  phenomena  included  in  the 
model.  The  questions  raised  in  the  Experimental  Approach  section  of 
this  report  will  be  answered  after  comparing  the  fatigue  test  results 
with  the  model  output.  Completion  of  the  research  objectives  are 
dependent  on  these  answers. 


DYNAMIC  MATHEMATICAL  MODEL 

The  purpose  of  the  dynamic  mathematical  model  is  to  answer  the 
questions  posed  in  the  Experimental  Approach  section  of  this  report. 

The  model  will  be  used  to  validate  the  technical  approach  through 
analysis  of  the  1018  and  4340  steel  fatigue  test  results.  Analysis  of 
the  4340  steel  results  will  explain,  through  the  model,  the  steady 
increase  in  damping  that  cannot  be  represented  by  dislocation  phenomena. 

Beam  Models 

The  two  dynamic  mathematical  models  developed  to  analyze  damping 
are  based  on  the  Euler-Bernoul 1 i  beam  equation  (Ref  72).  These  models 
are  viscously  damped  with  the  damping  based  on  a  dislocation  damping 
expression,  including  both  strain/ampl itude-dependent  and  frequency- 
dependent  damping.  In  addition,  the  fatigue  specimen  model  simulates 
crack  growth  with  a  torsional  spring,  kc(a),  that  is  a  function  of  crack 
length,  a. 

Figures  24  and  25  illustrate  models  of  the  four-point  beam  fatigue 
and  symmetric  cantilever  beam  damping  specimens,  respectively.  In  both 
figures,  m  is  the  mass  of  the  beam  per  unit  length,  E  is  Young's  modulus, 
h  is  the  beam  height,  b  is  the  beam  thickness,  1  is  the  beam  length,  and 
x  and  z  are  the  beam  coordinates.  In  Figure  24,  k^e  is  the  torsional 
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stiffness  of  the  knife-edge  support,  ^(a)  is  the  torsional  stiffness  of 
the  fatigue  crack  (as  described  above),  m  is  the  mass  of  the  exciter, 

1  is  the  distance  from  the  knife-edge  support  to  the  exciter,  and 

F  (x,t)  is  the  time-dependent  forcing  function  of  the  exciter.  In 
Figure  25,  w^(t)  is  the  time-dependent  acceleration  of  the  base  support 
for  the  symmetric  cantilever  beam. 

Symmetric  Cantilever  Beam  Damping  Specimen.  Appendix  F  derives  the 
equation  of  motion  and  mass,  stiffness,  and  damping  terms  for  the  symmetric 
cantilever  beam  model.  Results  from  the  internal  damping  tests  are  analyzed 
using  Equation  F-13  from  Appendix  F.  This  equation  determines  the  damping 
ratio,  at  resonant  frequency,  from  the  beam  acceleration  and  strain. 

Four-Point  Beam  Fatigue  Specimen.  Appendix  G  derives  the  boundary 
conditions,  equation  of  motion  and  mass,  stiffness,  and  damping  terms 
for  the  four-point  beam  fatigue  model.  Damping  data  calculations  for 
the  fatigue  tests  are  programmed  into  the  test  control  computer.  These 
calculations  are  based  upon  Equations  F-13  and  G-8  in  the  appendixes. 

The  following  section  develops  an  expression  for  the  frequency- 
dependent  decrement  that  must  be  determined  to  isolate  the  fatigue- 
induced  damping. 

Frequency-Dependent  Decrement 

A  simple  expression  for  frequency-dependent  decrement,  A^,  is 
provided  by  Zener  (Ref  1,  73)  for  a  beam  vibrating  in  simple  harmonic 
flexure  at  a  frequency  w.  According  to  the  Zener  thermoelastic  theory, 
the.  frequency-dependent  decrement  for  a  beam  is: 
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where  t,  the  relaxation  time  for  heat  flow  across  a  rectangular  beam,  is 
given  by: 
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and  a  is  the  linear  expansion  coefficient,  T  is  absolute  temperature,  c 
is  the  specific  heat  per  unit  volume,  and  k  is  the  thermal  conductivity. 
The  decrement  given  by  Equation  40  is  maximum  at  a  frequency  w  =  1/x  and 
falls  off  gradually  to  zero  for  very  high  and  very  low  frequencies. 

Based  on  data  presented  in  Table  4  (Ref  74),  Figure  26  reflects  the 
relaxation  damping  decrement  of  a  beam  modeled  by  Equation  40.  The 
curve  of  Figure  26  describes  a  single  relaxation  mechanism.  There  are  a 
number  of  such  mechanisms,  each  with  its  own  relaxation  time,  x.  Through 
linearity,  these  mechanisms  can  be  superimposed  to  provide  the  total 
damping  decrement  (Ref  71). 

Comparison  of  Fatigue  Test  and  Mathematical  Model 

In  summary,  the  steps  used  to  compute  the  estimated  internal 
damping  for  a  four-point  beam  model  undergoing  fatigue  crack  growth  are 
presented  in  Table  5.  In  applying  these  steps  to  1018  steel,  the  fol¬ 
lowing  assumptions  were  made  concerning  material  constants  (refer  to 
Table  6  for  a  listing  of  material-dependent  constants):  (1)  use  data 
previously  collected  to  estimate  constants,  (2)  if  immediate  data  are 
unavailable,  use  data  from  literature,  and  (3)  if  literature  data  are 
unavailable,  estimate  from  similar  materials.  Refer  to  Appendixes  C  and 
D  for  an  application  of  the  steps  presented  in  Table  5. 

True  stress-strain  for  1018  steel  is  modeled  by  Equation  34  (Figure 

3 

27)  with  k  =  84.2  x  10  psi  and  n  =  0.110.  Figure  28  is  a  plot  of  the 
P 

output  from  a  computer  model  depicting  a  1018  steel  beam  subjected  to 
fatigue  loading.  The  model,  outlined  in  Appendix  C  with  data  from  Table 
6,  matches  the  damping  characteristics  of  the  fatigue  specimen  as  the 
crack  grows.  The  sharp  increase  in  damping,  demonstrated  by  the  speci¬ 
men  for  a  large  plastic  zone,  is  also  exhibited  in  the  computer  generated 
model . 

Figure  29  is  a  plot  of  the  results  of  a  computer  model  depicting  a 
4340  steel  beam  subjected  to  fatigue  loading.  Data  for  the  computer  are 
shown  in  Table  6.  This  plot  is  compared  with  Figure  16.  The  model  has 


a  nonzero  value  for  the  constant  c  This  damping  effect  allows  the 
model  to  accurately  reflect  the  damping  conditions  of  the  steel  specimen. 
Differences  are  observed  at  12,000  and  23,000  cycles  where  sharp  damping 
peaks  exist  in  the  beam  specimen  output.  An  initial  impression  of  these 
peaks  is  that  they  are  due  to  external  causes,  as  they  appear  in  all  the 
4340  steel  specimens  at  consistent  frequencies. 

The  mathematical  model  has  many  simplifications.  The  plastic  zone 
is  treated  as  a  square  area,  where  each  side  is  twice  the  plastic  zone 
radius  in  length.  This  is  a  simple  treatment  with  room  for  refinement. 
Strain  at  the  crack  tip,  during  damping  data  collection  (Equation  13), 
is  assumed  to  be  a  function  of  the  reduced  cross  section  at  the  crack. 
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This  assumption  discounts  a  large  stress  concentration  that  occurs  at 
the  crack  tip  for  elastic  behavior  (Ref  60).  The  disregard  for  stress 
concentration  is  deemed  acceptable  in  this  instance  because  the  plastic 
yielding  that  occurs  within  the  plastic  zone  greatly  mitigates  stress 
concentration.  To  completely  ignore  the  stress  concentration  is  not 
realistic,  but  justifiable  for  an  initial  modeling  attempt.  The  ratio 
of  maximum  stress  to  nominal  stress  is  less  than  the  elastic  stress  con¬ 
centration  factor,  and  can  be  interpreted  as  an  elastic-plastic  stress 
concentration  factor,  which  must  be  determined  experimentally. 


CONCLUSIONS 

Conclusions  about  the  main  objective,  that  is,  the  crack 
growth/ internal  damping  relationship,  and  conclusions  about  the 
mathematical  model  used  in  the  analysis  of  the  damping  and  fatigue 
test  results  are  presented. 
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Damping  and  Fatigue  Tests 


The  major  objective  to  develop  a  relationship  between  crack  propaga¬ 
tion  and  internal  damping  was  accomplished.  Fatigue  tests  were  conducted 
for  validation  of  the  relationship.  Conclusions  based  on  the  test  results 
are  presented  and  the  reliability  of  the  damping  measurements  made  during 
the  fatigue  tests  is  addressed  in  the  following. 

Damping  Mechanisms.  Observations  of  fatigue  specimens  isolated  two 
distinct  damping  mechanisms.  The  first  was  classified  as  a  dislocation 
damping  mechanism  and  was  observable  in  1018  steel  specimens  only  after 
extensive  fatigue  crack  growth  under  constant  load  produced  a  large  plas¬ 
tic  region.  Dislocation  damping  is  well  documented  (Ref  33  through  43), 
and  the  results  of  the  1018  steel  beam  fatigue  tests  reflect  that  mecha¬ 
nism.  Observation  of  this  mechanism  was  an  objective  of  the  4340  steel 
beam  fatigue  tests.  However,  for  the  4340  specimens,  fatigued  under  a 
constant  stress  intensity  factor,  the  plastic  zone  radius  was  too  small 
(on  the  order  of  0.003  inch)  to  furnish  a  measurable  damping  increase. 
Indeed,  any  damping  increase  generated  by  a  plastic  zone  was  obscured  by 
a  second  damping  mechanism  that  was  observable  within  a  few  cycles  of 
crack  initiation.  The  damping  ratio  produced  by  this  mechanism  is 
directly  proportional  to  the  maximum  beam  rotation  at  the  crack  loca¬ 
tion.  The  second  damping  mechanism  was  classified  as  viscous  damping 
expressed  by  Equation  39  and  was  attributed  to  the  rubbing  together  of 
the  two  crack  surfaces  as  the  specimen  vibrated. 

No  explanation  is  immediately  available  for  the  cause  of  the  sharp 
peaks  appearing  in  the  damping  plots  of  the  4340  steel  beam  specimens. 
Nevertheless,  the  occurrence  of  the  peaks  at  consistent  frequencies 
suggests  the  source  is  external  to  the  specimens.  This  conclusion  is 
based  on  the  knowledge  that  each  specimen  had  a  different  initial 
resonant  frequency.  Knife-edge  stiffness  is  calculated  from  resonant 
frequency,  thus  all  the  4340  beam  specimens  had  different  knife-edge 
stiffnesses.  The  crack  length  and  general  beam  mode  shapes  are  thus 
different  for  each  beam  at  the  frequencies  where  the  damping  peak 
occurs . 


43 


Test  Procedures.  Results  from  the  symmetric  beam  damping  tests 
indicate  that  the  strain  gage  lead  wire  contribution  to  damping  is  stable 
and  produces  repeatable  results  that  may  be  treated  as  an  external  damping 
source.  For  this  reason,  lead  wire  damping  is  assumed  to  be  part  of  the 
damping  datum  established  for  each  fatigue  test. 


Mathematical  Model 

The  mathematical  computer  model  developed  to  express  damping  in 
materials  subject  to  crack  growth  is  a  combination  of  analytical  and 
empirical  relations.  The  empirical  relations  are  primarily  associated 
with  the  viscous  damping  observed  in  the  4340  steel  beam  specimens. 

Deduced  as  friction  between  the  two  crack  surfaces,  this  viscous  damping 
is  repeatable  in  4340  steel  beam  specimens.  Another  empirical  relation, 
expressed  in  F.quation  30,  is  assumed  to  exist  between  the  torsional  stiff¬ 
ness  at  the  crack  and  the  crack  length.  This  relationship  closely  resem¬ 
bles  a  moment  of  inertia  expression  and  gave  consistent  results  for  all 
four  4340  beam  specimens. 

The  analytical  expressions  used  to  model  dislocation  damping  and 
plastic  strain,  that  is,  F,quations  13  and  21,  accurately  modeled  large 
damping  changes  in  the  1018  steel  specimens.  This  conclusion  reflects 
the  difficulty  in  accurately  measuring  small  changes  in  damping  caused 
by  dislocation  damping,  while  acknowledging  the  model's  credibility  for 
simulating  large  damping  changes  which  occurred  after  many  thousands  of 
fatigue  cycles. 


RECOMMENDATIONS 

Additional  testing  is  required  to  confirm  the  mathematical  model 
and  the  viscous  damping  observed  in  the  4340  beam  specimen.  Other 
materials  must  be  fatigued  using  both  constant  load  and  stress  intensity 
factor.  Greater  accuracy  is  required  in  the  measurement  of  damping. 
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Damping  ratios  in  the  4340  steel  beam  specimens  could  be  measured  to  50 
x  10  Damping  changes  on  the  order  of  3  x  10  ^  are  suggested  by  the 
mathematical  model  for  simulating  dislocation  multiplication  in  the 
plastic  zone. 

The  mathematical  model  should  be  refined  to  better  represent  the 
plastic  zone  in  a  beam.  The  plastic  zone  shape,  as  handled  by  the  model, 
is  shown  in  Figure  30.  Using  a  square  shape  for  the  plastic  zone  is  a 
gross  approximation,  particularly  for  the  1018  steel  specimens.  Under 
fatigue  loading,  the  1018  steel  specimens  displayed  plastic  hinges  at 
the  notch  location.  Figure  31  illustrates  a  plastic  hinge,  as  observed 
in  the  1018  steel  specimens.  The  simple  simulation  used  in  the  mathe¬ 
matical  model  does  not  accurately  depict  the  plastic  material  opposite 
the  cracked  surface.  A  refinement  that  accounts  for  the  plastic  hinge 
effect  would  directly  influence  the  calculated  damping  with  the  plastic 
zone  and  improve  the  model's  accuracy.  Along  with  shape  is  the  question 
of  strain  within  the  plastic  zone  during  damping  data  collection.  Dis¬ 
location  damping  is  a  function  of  strain,  and  when  the  beam  is  vibrated 
to  collect  damping  data,  strain  within  the  plastic  zone  is  affected  by 
both  a  reduced  cross  section  and  plastic  flow.  The  model  ignores  stress 
concentrations  at  the  crack  tip  on  the  basis  that  the  concentration  is 
greatly  diminished  by  the  presence  of  a  plastic  region  (Ref  11).  Further 
work  is  required  in  this  area  to  ascertain  what  effect  each  phenomenon 
has  upon  the  strain  within  the  plastic  zone. 

Experimental  data  should  also  be  collected  to  determine  the  elastic- 
plastic  strain  distribution  around  the  notch  tip.  Equation  41  could 
then  be  refined  to  include  an  elastic-plastic  stress  concentration  fac¬ 
tor,  reflecting  plastic  flow  at  the  tip  and  accurately  estimating  the 
plastic  strain  in  the  plastic  zone. 

Damping  has  potential  as  a  crack  growth  monitoring  tool  where 
external  damping  can  be  tightly  controlled.  Isolation  from  external 
effects  is  possible  in  a  laboratory  where  conditions  are  rigidly  con¬ 
trolled.  Outside  a  laboratory,  changes  in  external  damping  will  often 
overwhelm  smaller  crack  growth  induced  damping.  A  technique  is  needed 
to  localize  or  "zoom-in"  on  crack- induced  damping.  Examination  of 
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localized  displacements  (on  the  order  of  0.5  microns)  is  possible  with 
the  Sharpe  technique  (Ref  75).  This  technique  also  has  the  advantage  of 
eliminating  strain  gages  and  their  affiliated  wires  and  damping  affects. 
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LIST  OF  SYMBOLS 


6(x-$) 


Linear  expansion  coefficient  and  constant 

Material-dependent,  crack  growth  rate  constant 

Unit  impulse  or  delta  function  at  the  location  £ 

Total  decrement,  Au  +  A 
n  I 

Average  specific  damping  capacity 

— c  — c 

Average  corrected  total  decrement,  A  +  A 

H  I 

External  damping 
Amplitude-dependent  decrement 

Povolo's  corrected  amplitude -dependent  decrement 
Average  amplitude-dependent  decrement 

Average  Povolo's  corrected  amplitude-dependent  decrement 
Frequency-dependent  decrement 

Average  corrected  frequency-dependent  decrement 
Frequency-dependent  decrement  for  the  i^  mode 


Difference  between  the  maximum  and  minimum  values  of  the 
stress  intensity  factor  during  fatigue  cycle,  Kmflx  * 

Internal  damping  due  to  plastic  strain 

4(1-v)/tt2 


Energy  dissipated  per  cycle  in  a  volume  element 


Damping  decrement  for  beam  specimen  with  no  sets  of 
wires  removed 


Damping  decrement  with  one  set  of  strain  gage  wires  removed 
Damping  decrement  with  two  sets  of  strain  gage  wires  removed 
Damping  decrement  with  three  sets  of  wires  removed 
Damping  decrement  with  all  strain  gage  wires  removed 
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E(Xo) 


Damping  decrement  of  i  mode  with  all  wires  removed 
Crack  growth  during  fatigue  cycle 
Cottrell  misfit  factor 

Fatigue  cycle  elastic  strain  at  maximum  load 

Fatigue  cycle  strain  at  minimum  load 

Plastic  strain 

Plastic  strain  per  cycle 

Plastic  strain  for  fatigue  cycle 

Static  plastic  strain  fracture  value  at  1/2  cycle 

Ultimate  strain  for  static  fracture 

Critical  strain  above  which  plastic  strain  occurs 

Average  value  of  the  cyclic  plastic-strain  range 

ln(A  /A.),  true  strain 
o  1 

Strain  in  beam  at  location  x 

o 

Cyclic  strain  amplitude 

Viscous  damping  ratio  for  the  i^  mode,  =  C^/C^ 

Phase  angle,  8  -*■  it/ 2  as  w  -*  u . 

A"/A 

Eigenvalue  of  i  ^  mode  of  vibrating  beam 

Total  length  of  movable  dislocation  line  per  unit  volume 
(dislocation  density) 

Dislocation  density  in  elastic  region 

Immobile  dislocation  density 

Dislocation  density  in  plastic  zone 

Poisson's  ratio 

Location  of  unit  impulse  function 
Dens ity 
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0 


Nominal  outer  fiber  stress 


0 

a 


a . 

10 


o 


t 


a 

u 

o(0) 


i 

* 


<^i(x) 
1>i  '  (x) 

*."(X) 


X 


4> 

w 


w . 
1 


w 


n  1 


Q 


a 


a 


0 


a 


j-1 


a  . 
J 

A 


A. 


A 

o 


Maximum  cyclic  stress 
Stress  constant 
P/A^,  true  strain 

Ultimate  stress  for  static  fracture 
Stress  at  zero  plastic  strain 

2  2 

Relaxation  time  for  heat  flow  across  a  beam,  t  =  h  c/ir  k 

Orientation  factor  (0.5  and  0.32  for  b.c.c.  and  f.c.c., 
respectively) 

Characteristic  shape  function  associated  with  i^  mode 
of  vibrating  beam 

Derivative  of  i*"*1  characteristic  shape  function  with 
respect  to  x,  describes  slope  of  beam 

Second  derivative  of  i^  characteristic  shape  function 
with  respect  to  x,  describes  curvature  of  beam 

Strain  independent  probability  for  the  remobilization  of 
immobile  dislocations 

Strength  parameter 

Frequency  of  forcing  function 

Resonant  frequency  of  the  it*1  mode,  =  (K^/M^)*^ 

Fundamental  beam  frequency 

Orientation  factor 

Crack  length 

Initial  crack  length 

Crack  length  after  i^-l  fatigue  cycle 
Crack  length  after  fatigue  cycle. 

2 

Dimensionless  constant  and  for  the  G-L  Theory  A  =  trpb 
Instantaneous  cross-sectional  area  under  a  given  load 
Original  cross-sect ional  area  before  any  load  is  applied 
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A’ 


l/$bKf 


Material  constant 

b  Burger's  vector  and  beam  thickness 

B  Material  temperature-dependent  term 

B.  Amplitude  scaling  factor  associated  with  i^  mode  of 

vibrating  beam 

b  Material  constant  relating  crack  length  to  beam 

torsional  stiffness  at  crack 

2 

C  Material  constant;  for  the  G-L  Theory  C  =  2Gb  /ir(l-v); 

and  for  the  Manson-Cof f in  law,  the  ductility  exponent 

C1  C40LN3/’Lc 

Cn  Kefa/L 

z  c 

J2A  I,  4Bw/tt2C 
3  0  c 

Generalized  damping  of  i^  mode 
C  ^  Generalized  critical  damping  of  i^  mode 

c  Specific  heat  per  unit  volume 

c  ^  Crack  damping  coefficient 

c  Equivalent  viscous  damping 

c(x)  Beam  damping  coefficient 

1/2 

d  Dislocation  cell  size,  d  =  K/A 

dw  Plastic  energy  absorbed  per  cycle 

D %  Fatigue  life  damage  assessment  in  percent,  100%  (N^/N  ) 

F.  Young's  modulus  of  beam  material 

Ej  Exponential  integral  function  of  first  order 

Exponential  integral  function 
f  Constant  of  order  1 

F(a/b)  1.122  -  1 . 4( a/b)  +  7.33(a/b)2  -  13.08(a/b)3  +  lA.0(a/b)4 
f  Amplitude  scaling  factor  of  exciter  forcing  function 
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F  (x,t) 
ex 

G 

h 


I 

j 

k 

K 


ave 

kc(a) 


Ic 


ke 


max 


min 


1 

I. 

1 


e 

1 

ex 

I,F 

I*., 


Exciter  forcing  function 
Shear  modulus 
Height  of  beam 
Mode  number 

3 

Moment  of  inertia  of  beam,  I  =  th  / 1 2 
Fatigue  cycle  number 

Thermal  conductivity;  for  Manson-Cof f in  law,  1/C 

Orientation  and  anisotropy  parameter  and  material 
constant,  K  =  20  for  Fe 

Average  stress  intensity  factor  over  fatigue  test 

Torsional  spring  stiffness  of  crack 

Generalized  stiffness  of  i^k  mode 

Stress  intensity  factor 

Critical  stress  intensity  factor 

Torsional  spring  stiffness  of  knife-edge  support 

Maximum  value  of  stress  intensity  factor  during  fatigue  cycle 

Minimum  value  of  stress  intensity  factor  during  fatigue  cycle 

Strength  coefficient 

Length  of  beam 

Mobile  dislocation  density 

Fatigue  specimen  distance  from  knife  edge  to  loading  point 

Impurities  loop  length 

Effective  length  of  dislocation  loop 

Distance  from  knife-edge  support  to  exciter 

Fraction  of  fatigue  life  expended 

Network  loop  length 


m 


Hysteresis  energy  exponent 


56 


M  Moment  on  beam  section 

M  Moment  due  to  crack 

c 

M^  Generalized  mass  of  i^*1  mode 

M^^  Moment  due  to  knife-edge  support 

m(x)  Mass  per  unit  length  of  beam 


M(x) 

n 


N 

N 


N, 


N, 

fp 

N* 


Moment  along  beam,  M(x)  =  El  $  "(x) 

Strain  hardening  exponent 

Fatigue  cycle  number 

Number  of  fatigue  cycles  completed 

Number  of  cycles  to  complete  fracture 

Predicted  number  of  cycles  to  complete  fracture 

Number  of  cycles  in  a  stage  of  fatigue  life 


P 


R, 

i 

r 

P 

s 


SD 


t 

T 

TF 

U 


Fatigue  or  static  load 
Generalized  coordinate  of  it^1  mode 

First  derivative  of  generalized  coordinate  with  respect 
to  time 

Second  derivative  of  generalized  coordinate  with  respect 
to  time 

Characteristic  shape  function  ratio  of  nonzero  constants 
Plastic  zone  radius 
Dislocation  mean  free  path 
Standard  deviation 

Time  variable  or  thickness  of  beam 
Absolute  temperature 

Transfer  function  amplitude  ratio  at  u  , 

n  1 

Integrated  plastic  work  expended  in  the  plastic  zone  or 
peak  potential  energy  per  cycle 


Uj  Rate  of  immobilization  of  mobi  l^dis  locat  i  ons  or  peak 

potential  energy  per  cycle  of  i  mode 


>;  ' 


V 


Specimen  volume 


Elastic  region  volume 
Plastic  zone  volume 


V 

e 

V 

P 

W  Amplitude  scaling  factor  for  beam  support  displacement 

w^(t)  Time-dependent  displacement  of  beam  support 

w^(t)  Time-dependent  acceleration  of  beam  support 

Energy  dissipated  per  cycle 

Energy  dissipated  per  cycle  of  i^  mode 

W^.  Total  hysteresis  energy  to  cyclic  failure 

Total  hysteresis  energy  to  static  failure 

w(x,t)  Time-  and  location-dependent  displacement  of  beam 

w(x,t)  Time-  and  location-dependent  velocity  of  beam 

w(x,t)  Time-  and  location-dependent  acceleration  of  beam 

w'(x,t)  Derivative  of  w(x,t)  with  respect  to  x 

w"(x,t)  Second  derivative  of  w(x,t)  with  respect  to  x 

w'"(x,t)  Third  derivative  of  w(x,t)  with  respect  to  x 

w  Peak  potential  energy  per  cycle  in  a  volume  element 

x  Coordinate  of  longitudinal  axis  of  beam 

x  Known  location  on  beam 

o 

z^(t)  Forcing  function  for  generalized  coordinate  equation  of 

motion 

z^(t)  Displacement  function  for  i^  generalized  coordinate 
equation  of  motion 

Zj,  Aigjalitude  scaling  factor  for  displacement  function  of 

i  mode 


Table  1.  Characteristics  of  SAE  4340  Steel 


Item 

Description 

Material 

SAE  4340  Steel 

Composition 

C  -  0.40,  Mn  -  0.81,  P  -  0.018, 

S  -  0.019,  Si  -  0.25,  Ni  -  1.73, 

Cr  -  0.87,  Mo  -  0.24 

Heat  Treatment 

The  fatigue,  tension,  compact  tension  and 
damping  specimens  were  austenitized  at 
1525°F  in  a  neutral  salt  bath,  quenched 
in  still  oil  and  room  temperature  and 
tempered  at  1100°F  for  2  hours. 

Engineering 

Tensile  Properties 

Average  tensile  properties  for  three 
specimens  are  given  below: 

Engineering  0.2%  offset, 

yield  strength  (psi)  -  141,300 

Engineering  ultimate 

strength  (psi)  -  150,200 

Hardness  Rockwell  C  -  37 

Table  2.  Summary  of  Static  Test.  Results 


Parameters 

Units 

Values 

. - - -  ‘ 

-  '  '  -  '  “ - - - - 

•  ---  - 

s 

u 

ps  i 

225,000 

t 

pu 

in .  /  in . 

0.748 

W 

u 

lb-in. /in. ^  (Equation  37) 

153,800 

w 

3 

lb- in. /in.  ( Equation  38 r 

154,000 

U 

.  .  '  . 
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Table  3.  Fatigue  Specimens  at  Failure 


Table  4.  Values  Used  in  Equation  40  to  Obtain 
the  Curve  Shown  in  Figure  26 
(Ref  74)  for  Steel 


Table  5.  Computation  of  Internal  Damping  in  Fatigue  Model 


Step 


Drescription 


1 

2 

2.1 

2.2 


2.3 


2.4 


2.5 


2.6 


2.7 

2.8 
2.9 

2.  10 


Define  material  and  crack  growth  constants 

Repeat  step  2  for  fatigue  cycle  j  =  1  to  N 

Using  the  crack  growth  constant,  o  ,  and  crack  length  from 
the  previous  fatigue  cycle,  a  ,  compute  the  crack  length 
for  the  current  fatigue  cycleJa^  (Equations  31  and  G-10) 

From  known  values  of  the  static  plastic  strain  fracture 
value,  t  u>  the  final  crack  length,  a^,  and  the  ductility 
exponentpUC,  compute  the  estimeted  plastic  strain,  z  .  in 
the  current  fatigue  cycle  (Equation  30) 

Compute  the  fraction  of  fatigue  life  expended,  LF 
(Equation  29) 

Determine  the  plastic  zone  radius,  r  ,  for  the  current 
fatigue  cycle  (Equation  11  or  12)  p~* 

Calculate  the  crack  torsional  stiffness,  k  .,  for  the 


current  fatigue  cycle  (Equation  G-12) 


cj 


With  the  knife-edge  stiffness,  k^  ,  determine  the 
eigenvalue,  eigenvector,  $j(x^,  and  resonant  fre¬ 

quency,  w^,  or  the  fundamental  mode  of  vibration 
(Equation  G-4) 

Determine  amplitude  scaling  factor,  (Equation  F-8) 

Calculate  plastic  zone  dislocation  density  (Equation  21) 

— c 

Compute  the  average  amplitude-dependent  damping,  A  ,  from 
the  vibration  strain  level,  e^  (Equation  13)  J 

Determine  the  frequency-dependent  damping,  Aj  (Equation 
F-16) 

stop. 
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Table  6.  Material -Dependent  Constants  for  1018  and  4340  Steel 


Constant 

Description 

Units 

1018 

4340 

Ai 

Dislocation  density  to  plastic 
strain  coefficient 

disl/in. 2 

2.41  x  106 

6. 13  x  106 

b 

c 

Crack  length  to  beam  torsional 
stiffness  coefficient 

lb/ in. 

124  x  103 

124  x  103 

c 

Manson-Cof f in  ductility  exponent 

646  x  10‘3 

690  x  10‘3 

ci 

Internal  damping  coefficient 

50  x  10‘3 

50  x  10“3 

C2 

Internal  damping  coefficient 

in . / in . 

150  x  10'6 

100  x  10‘6 

c 

Specific  heat  per  unit  volume 

ps i/°F 

245 

245 

El 

Young's  modulus  x  beam  moment 
of  inertia 

lb/in.2 

300  x  103 

300  x  103 

k 

Thermal  conductivity 

6.5 

6.5 

a 

Expansion  coefficient 

in. /in. °F 

6.7  x  10"6 

6.7  x  10'6 

a 

c 

Crack  growth  rate  coefficient 

in.  /lb 
cycle 

16  x  10"24 

9.66  x  10'24 

e 

pu 

Static  plastic  fracture  value 
at  1/2  cycle 

in .  /  in . 

550  x  10*3 

748  x  10‘3 

Ao 

Dislocation  density  for  plastic 
strain  equal  to  zero 

disl/ in. 2 

1  x  108 

1  x  108 

a 

y 

Yield  stress 

psi 

3 

62  x  10 

144  x  103 

X 

Strain  independent  probability 
for  the  remobilization  of 
immobile  dislocations 

10.7 

8.61 

62 
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Figure  1.  Edge  dislocation. 
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Figure  3 


Formation  of  a  microcrack  by  a  combination  of  (three) 
dislocations  (from  Yokobori,  Ref  11). 


65 


Strain  Amplitude  X  Damping  Ratio  X  10 


(Strain  Amplitude)  1  X  10  3 


Figure  4.  Granato-Lucke  plot  of  a  1018  steel  cantilever  beam. 
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Figure  5.  Logarithmic  stress-plastic  strain  curve  (from  Feltner 
and  Morrow,  Ref  47). 
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stress-strain  curve:  4340  steel  (from  Feltner 
lorrow,  Ref  47). 
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Average  Plast 


applied  load 


CL 


Figure 


8.  Successive  stages  of  ductile 
failure  (from  Blauel,  Ref  12) 


Figure  9.  Fatigue  test  equipment  schematic. 
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Figure  10.  Fatigue  test  beam  specimen. 
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Figure  11.  Damping  test  equipment  schematic. 


Figure  12.  Damping  test  beam  specimen. 
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Figure  13.  True  stress-strain  curve:  4340  steel. 
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Figure  14.  Strain  amplitude  dependent  dislocation  damping. 
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Figure  15.  Beam  Specimen  A:  Crack  extension,  resonant 
frequency  and  damping  ratio. 
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Figure  16.  Beam  Specimen  B:  Crack  extension,  resonant 
frequency  and  damping  ratio. 
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Figure  18.  Beam  Specimen  D:  Crack  extension,  resonant 
frequency  and  damping  ratio. 
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Figure  19.  Beam  Specimen  B:  Damping  versus  beam  slope. 
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Figure  20.  Beam  Specimen  C:  Damping  versus  beam  slope. 


Beam  Specimen  D:  Damping  versus  beam  slope. 
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1018  STEEL  BEAM  SPECIMEN  (SMOOTH) 
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23.  1018  steel  beam  specimen:  Crack  extension,  resonant 

frequency  and  damping  ratio. 
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Figure  26.  Frequency  dependent  damping,  experimental  results. 
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stress-strain  curve:  1018  steel 


COMPUTER  SIMULATION  OP  CRACK  GROWTH  IN  1018  STL 
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Figure  28.  1018  steel  beam  simulation:  Crack  extension,  resonant 

frequency  and  damping  ratio. 


COMPUTER  SIMULATION  OF  CRACK  GROWTH  IN  SEAM  B 
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frequency  and  damping  ratio. 


Beam  Notch 


Figure  31.  Plastic  hinge 


Appendix  A 


FATIGUE  TEST  PROCEDURE 


This  appendix  presents  a  detailed  description  of  the  procedure  used 
to  fatigue  test  1018  and  4340  steel  beam  specimens. 


TEST  EQUIPMENT 

Table  A-l  contains  a  list  of  equipment  used  in  the  fatigue  test. 

Refer  to  Figure  9  in  the  main  text  of  this  report  for  a  general  arrange¬ 
ment  of  the  equipment  listed  in  Table  A-l. 

A  vibration  generator  excites  a  fatigue  specimen  to  collect  strain 
and  force  data  used  in  the  calculation  of  internal  damping  for  the  beam. 

An  impedance  head  located  between  the  vibration  generator  and  specimen 
provides  the  force  data  for  the  damping  calculation.  Three  strain  gages 
are  attached  to  the  specimen  (see  Figure  10).  These  gages  provide  strain 
measurements  during  fatigue  loading  and  vibration  generator  excitation. 

The  specimen  is  vibrated  at  a  frequency  selected  via  computer  control  on 
the  frequency  synthesizer.  The  frequency  synthesizer  output  signal  is 
boosted  by  a  power  amplifier  to  provide  strains  of  measurable  amplitude. 

The  Lab-Datax  Computer  controls  the  fatigue  load  of  the  MTS  test 
machine  and  performs  the  calculations  to  determine  internal  damping. 

Crack  propagation  gages  connected  to  an  analog  to  digital  port  on  the 
computer  monitor  the  fatigue  crack  length  that  the  computer  uses  to 
maintain  a  constant  stress  intensity  factor.  The  computer  also  controls 
the  number  of  fatigue  cycles  between  damping  data  cycles  and  the  frequency 
range  selected  on  the  frequency  synthesizer. 


A-l 


Table  A-l.  Fatigue  Test  Equipment 


1  -  Data  Translation  Lab-Datax  Computer  w/Hard  Disk  Drive,  D/A 
Converter,  and  A/D  Converter 

1  -  Scientific  Atlanta  Model  SD375  Spectrum  Analyzer 

1  -  Wilcoxon  Research  Model  F5B  Vibration  Generator 

1  -  Wilcoxon  Research  Model  Z12  Impedance  Head 

1  -  Accelerometer  -  Force  Amplifier  Unit 

1  -  Validyne  Model  MC1-3  Amplifier  Enclosure 

2  -  Validyne  Model  SG71  Strain  Gage  Amplifier 

1  -  Hewlett  Packard  Model  3325A  Frequency  Synthesizer 

1  -  Hewlett  Packard  Model  3495A  Scanner  w/Optional  001  Relay 

Assembly 

1  -  Hewlett  Packard  Model  467A  Power  Amplifier 

1  -  MTS  Test  Machine 

2  -  Micro-Measurements  Model  TK-09-CPA0 1-005  Crack  Propagation 

Gage 

3  -  Micro-Measurements  Model  EA-06- 125BZ-350  Strain  Gage 


The  spectrum  analyzer  collects  strain  and  force  data  and  converts 
them  to  a  transfer  function  ratio.  The  transfer  function  is  transferred 
to  the  computer  where  the  damping  calculation  is  made  and  the  frequency 
range  selection  is  updated.  Strain  input  to  the  spectrum  analyzer  is 
under  computer  control  with  the  switching  function  handled  by  the  scanner. 


SPECIMENS 

The  1018  or  4340  steel  beam  specimen,  shown  in  Figure  10,  is  mounted 
in  the  MTS  test  machine  using  a  knife-edge  support  fixture.  The  specimen 
is  initially  uncracked,  with  a  small  notch  at  its  midpoint.  A  four-point 
load  is  applied  to  the  beam  by  the  hydraulic  ram  of  the  MTS  test  machine. 


A-2 


The  initial  fatigue  load,  ranging  from  1,250  to  1,400  pounds  for  the 

four  specimens  tested,  produce  a  maximum  stress  in  the  notch  of  160  ksi 

1/2 

and  a  maximum  stress  intensity  factor  of  36  ksi-in.  '  . 

For  the  heat  treatment  used,  the  4340  steel  has  a  yield  strength 

1/2 

of  141  ksi  and  a  critical  stress  intensity  factor  of  77  ksi-in.  '  . 

With  the  specimen  fatigued  at  a  constant  stress  intensity  factor,  the 

-3 

maximum  plastic  radius  is  estimated  at  3.5  x  10  inches,  as  calculated 
from  Equation  12  in  the  main  text  of  this  report. 

The  1018  steel  was  untreated,  having  a  yield  strength  of  62  ksi. 

The  1018  steel  is  too  soft  to  obtain  a  valid  critical  stress  intensity 
factor  and  no  attempt  was  made  to  determine  one.  At  800  pounds,  the 

.3 

uncracked  specimen  has  a  maximum  plastic  radius  estimated  at  7.2  x  10 
inches,  as  calculated  from  Equation  12. 


TEST  PROCEDURE 

After  mounting  the  specimen  in  the  fixture  on  the  MTS  test  machine, 
the  Lab-Datax  computer  program  is  initialized.  The  initialization  pro¬ 
cess  sets  the  number  of  fatigue  cycles,  crack  gage  parameters , and  exciter 
frequency  range.  After  initializing  the  system,  2,000  load  cycles,  100 
pounds  each,  are  applied  to  the  beam  specimen.  This  series  of  load  cycles 
is  used  to  establish  a  damping  datum  for  the  uncracked  specimen. 

For  each  load  cycle,  maximum  strain  and  load  are  saved  on  disk  by 
the  computer.  After  each  set  of  forty  load  cycles,  the  computer  program 
executes  a  damping  data  cycle.  A  damping  data  cycle  begins  when  the 
frequency  synthesizer  is  set  to  sweep  a  range  of  frequencies  starting  10 
Hz  below  and  ending  10  Hz  above  the  fundamental  resonant  frequency  of 
the  specimen.  As  the  vibration  generator  sweeps  the  frequency  range, 
exciting  the  specimen,  strain  gage  and  force  transducer  data  are  trans¬ 
ferred  to  the  spectrum  analyzer.  These  data  are  transformed  into  a 
transfer  function  consisting  of  two  components:  real  and  imaginary. 

Data  from  the  real  component  of  the  strain- force  transfer  function  are 
interpreted  in  terms  of  the  phase  shift  between  the  input  forcing  func¬ 
tion  and  output  strain.  The  damping  ratio,  £,  may  be  calculated  from 
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the  real  component  after  determining  the  frequencies  of  minimum  and 
maximum  magnitudes,  f^  and  f^,  respectively,  using  the  following  for¬ 
mula, 

2  2 
f  -  f 
2  1 

C  =  —^-5 - S"  (A- 1 ) 

2(f2  +  fl ’ 

Equation  A-l  was  not  used  for  calculating  the  damping  ratio  within  the 

2  2 

computer  program  because  the  term  f^  -  f^  is  a  small  difference  between 
large  numbers  and  has  poor  accuracy. 

A  better  technique  for  measuring  damping,  during  steady-state  vibra¬ 
tion,  is  to  measure  the  maximum  amplitude  of  the  imaginary  component  of 
the  transfer  function,  TFr  .  Inverting  TFt  and  dividing  by  two 

is  the  equivalent  of  Equation  F-13  in  Appendix  F,  which  is  equal  to  the 

damping  ratio.  TF  occurs  .at  the  resonant  frequency  of  the  beam, 

1  ,  max 

fr-  Steady-state  vibration  was  selected  instead  of  using  a  ring-down 
technique  because  the  damping  ratio  changes  with  vibration  amplitude  as 
the  specimen  rings  down. 

A  correction  factor  must  be  included  in  the  TFt  calculation. 

I  ,max 

The  spectrum  analyzer  collects  data  while  the  synthesizer  sweeps  discrete 
frequencies,  resulting  in  an  averaged  transfer  function,  rather  than  the 
true  value  at  each  discrete  frequency.  The  correction  factor  is  the 
inverse  of  the  number  of  data  sets  collected  per  frequency  sweep. 

After  establishing  the  damping  datum,  the  computer  applies  the 
full  fatigue  load  to  the  specimen.  As  before,  the  fatigue  load  is 
cycled  forty  times  between  damping  data  collection  cycles.  During  the 
load  cycle,  crack  length  measurements  are  made,  and  the  fatigue  load  is 
adjusted  to  maintain  a  constant  stress  intensity  factor.  This  procedure 
continues  until  the  crack  passes  beyond  the  crack  propagation  gages.  At 
this  point,  feedback  from  the  crack  propagation  gages  indicates  a  constant 
crack  length  and  the  system  produces  a  constant  fatigue  load  through  the 
remainder  of  the  test. 
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NOISE 


The  principal  sources  of  noise  in  the  fatigue  and  damping  tests  are 
fluorescent  lights  and  electomagnetic  fields  from  vibration  generators. 

In  the  case  of  fluorescent  lights,  the  spectrum  analyzer  filters  out  the 
noise  occurring  at  60  and  120  Hz.  For  the  vibration  generators,  tests 
showed  there  were  electomagnetic  effects  only  when  the  strain  gage  wires 
were  in  direct  contact  with  the  generator.  During  the  fatigue  and  damping 
tests,  care  was  taken  to  assure  that  the  strain  gage  wires  never  came 
into  direct  contact  with  a  vibration  generator. 
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Appendix  B 


DAMPING  TEST  PROCEDURE 

This  appendix  describes  the  procedure  used  to  collect  damping  data 
from  4340  steel  symmetrical  cantilever  beam  specimens. 


TEST  EQUIPMENT 

Table  B-l  contains  a  list  of  the  equipment  used  to  conduct  the 
damping  tests.  Refer  to  Figure  11  in  the  main  text  of  this  report  for  a 
schematic  of  the  test  equipment. 


Table  B-l.  Damping  Test  Equipment 
1  -  Scientific  Atlanta  Model  SD375  Spectrum  Analyzer 

1  -  Validyne  Model  MC1-3  Amplifier  Enclosure 

2  -  Validyne  Model  SG71  Strain  Gage  Amplifier 

1  -  Hewlett  Packard  Model  3325A  Frequency  Synthesizer 
1  -  Hewlett  Packard  Model  467A  Power  Amplifier 
1  -  MB  Electronics  Model  PM-50  Shaketable 
4  -  Micro-Measurements  Model  EA-06- 125BZ-350  Strain  Gage 


The  frequency  synthesizer  is  set  to  operate  the  shaketable  at  the 
resonant  frequency  of  the  symmetric  cantilever  beam  damping  specimen. 
Output  from  a  strain  gage  and  accelerometer  are  transmitted  to  the 
spectrum  analyzer  where  the  data  are  transformed  into  a  transfer 
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function.  The  amplitude  of  the  transfer  function  is  read  from  the 
analyzer  and  converted  to  a  damping  ratio  using  Equation  F-13  in 
Appendix  F. 

Amplitude  of  vibration  is  adjusted  through  the  power  amplifier. 
At  each  setting,  the  maximum  cyclic  strain  is  read  from  the  spectrum 
analyzer. 


SPECIMEN 

Figure  12  in  the  main  text  illustrates  the  symmetric  cantilever 
beam  specimen  used  in  the  damping  tests.  The  beam  material  is  4340 
steel  that  has  been  heat  treated  to  150  ksi  yield  strength. 


TEST  PROCEDURE 

A  damping  specimen  is  mounted  on  the  shaketable.  The  amplitude  of 
vibration  is  initialized  to  a  low  level  and  the  frequency  synthesizer  is 
adjusted  to  the  fundamental  resonant  frequency  of  the  specimen,  f  . 

From  the  spectrum  analyzer,  the  amplitudes  of  maximum  cyclic  strain,  eq, 
and  acceleration,  a^,  are  read.  The  ratio  of  a^  to  eq  is  a  transfer 
function.  To  pinpoint  the  exact  resonant  frequency  and  its  associated 
damping  ratio,  strain  and  acceleration  data  are  collected  from  frequencies 
0.05  Hz  above  and  below  f  .  Interpolation  is  used  on  these  three  data 
points  to  compute  the  beam's  damping  ratio.  The  damping  interpolation 
process  is  based  on  fitting  the  three  data  points  to  a  curve  generated 
by  Equation  F-ll. 

The  above  procedure  is  repeated  for  the  second  and  third  modes  of 
the  specimen.  Having  completed  the  damping  data  collection  for  the  first 
three  modes  of  the  specimen,  at  the  selected  strain  levels,  the  amplitude 
of  vibration  is  increased  and  the  procedure  repeated. 
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Appendix  C 

ALGORITHM  FOR  MODELING  INTERNAL  DAMPING  FROM  CRACK  GROWTH 

The  internal  damping  of  a  beam  undergoing  fatigue  crack  growth  is 
calculated  based  on  an  estimation  of  the  plastic  zone  size  and  material 
dependent  constants  related  to  crack  growth  rate  and  internal  damping. 
The  following  algorithm  is  used  to  calculate  the  internal  damping  of  a 
four-point  beam  with  a  propagating  crack. 

STEP  1.  Input  material  and  crack  growth  constants 

a^  Initial  crack  length  (in.) 

2 

Dislocation  density  coefficient  (disl/in.  ) 
af  Final  crack  length  (in.) 

a  Linear  expansion  coefficient  (in./in./°F) 

b  Material  constant  relating  crack  length  to 

torsional  stiffness  (lb-in.) 

C  Manson-Cof f in  ductility  exponent 

Internal  damping  coefficient 

Internal  damping  coefficient  (in. /in.) 

c  Specific  heat  per  unit  volume  (psi/°F) 

X  Strain  independent  probability  for  the 

remobilization  of  immobile  dislocations 

Young's  modulus  (psi) 

Young's  modulus  x  moment  of  inertia  (lb-in.  ) 

Static  plastic  fracture  value  at  1/2  cycle 
(in. /in. ) 

Ep  Maximum  strain  in  beam  during  damping  data 

collection  (In. /in.) 

a  Crack  growth  rate  constant  based  on  dislocation 

C  7  4 

theory  of  crack  growth  (in.  /lb  -cycle) 


E 

El 

e 

pu 


h 


Beam  height  (in.) 

Thermal  conductivity  (lb/sec-°F) 


k 


kke 

Torsional  stiffness  at  knife-edge  support 
(lb-in.) 

1 

Location  of  exciter  relative  to  knife-edge 

ex 

support  (in.) 

Xa 

Location  of  load  relative  to  knife-edge 
support  (in.) 

Ao 

2 

Dislocation  density  for  e^  =  0  (disl/in.  ) 

m 

Mass  of  beam  per  unit  length  (slug/ in.) 

m 

ex 

Mass  of  exciter  and  counterweight  (slug) 

P 

Fatigue  load  applied  to  beam  (lb) 

0 

y 

Yield  stress  (psi) 

STEP  2. 

Initialize  variables 

LF 

=  o., 

life  expended  variable 

flold 

tl 

o 

crack  length 

j 

=  i, 

,  2 

fatigue  cycle  counter 

t  = 

h  c 

v 

it  k 

relaxation  time  for  heat  flow 

STEP  3. 

while 

<aold  <=  V 

STEP 

3.  1 

Calculate  crack  length  for  current  cycle  based 

fanew  da 

~4 -  =  “c 

aold  V(a) 

Begin  by  calculating  the  stress  intensity  factor 
using  an  empirical  formula 

f (a/h)  =  1.122  -  1 . 4( a/h)  +  7.33(a/h)2  -  13.08(a/h)3  +  14.0(a/h)4, 
or 
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fU/h)  » 


1  /2  4 
'  0.923  +  0.199(1  -  sin») 


where  4>  = 

2.  n 


k  =  b  (h/a)  torsional  stiffness  at  crack 
c  c 


Rotation  of  beam  at  knife-edge  support 


( 2EI  +  1  k  )  1  P 
_  _ a  c  a 

"  -EI  (kke  +  kc)  ’ 


Nominal  bending  moment  between  load  points 


M  =1/2  1  P  ■  k,  0, 
nom  a  ke 


o  =  - 

nom  ,2 
h 


,  nominal  bending  stress 


Kt  =  a  /ira  f(a/h),  stress  intensity  factor. 
I  nom 


STEP  3.1.1  Assume  a  value  for  a 


STEP  3.1.2  Use  Newton-Raphson  Iteration  for  the  i 
iteration  of  a 


f  new  da  .  ,  , 

|  — ; — —  =  a  ,  residual  amount 

„  4.  .  c 

“old  KI  (a) 


R'  =  — r - 7— ,  slope  of  function 

KT  (a  A) 

I  new 

a  i+1  =  a  i  -  R/R' ,  i+lstiteration  of  a 
new  new  new 

Step  3.1  requires  numerical  integration  and  Newton-Raphson 
iteration.  This  case  has  one  equation  and  one  unknown  and 
should  not  require  extensive  computing. 

STEP  3.2  Sum  the  new  value  of  life  expended 


Aa  =  a  -  a  , . , 
new  old 


incremental  crack  growth 
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plastic  strain 


£ 

LF  +=  — ,  specimen  is  loaded  on  upstroke  only 
pu 

STEP  3.3  Compute  new  plastic  zone  size 

.  ,Ma  >? 
o  =  f  _I_a«L  , 

n  6ir  \  o 

y 

STEP  3.4  Calculate  new  torsional  stiffness  at  crack 


k  =  b  (h/a)3. 
c  c 


STEP  3.5  Compute  beam  fundamental  mode  characteristics 

including  eigenvalue,  x;  eigenvector,  $(x);  and  resonant 
frequency,  w.  The  six  unknowns,  x;  four  from  0(x):  A,  B, 
C,  and  D;  and  w,  are  computed  from  the  four  boundary 
conditions : 


w(0)  =  0, 


deflection  at  knife  edge  is 
zero 


Elw"(0)  =  k^^w'CO),  moment  at  knife  edge  is  equal 
e  to  the  torsional  stiffness  x 

s  lope 


EIw"(l/2)  =  k  w'(l/2),  moment  at  crack  is  equal  to 

the  torsional  stiffness  x 
s  lope 

w '  ’ ' (  1  / 2 )  =  0,  shear  at  crack  is  zero. 


Frequency  equation  (ignoring  the  small  amount  of  damping): 


w"Mi  -  K£  =  0, 


where 


=  f  EI[*  " 
In  i 


(x)j  dx, 


and 


M, 


=  fo*t2w "<*)  <*« +  *i2(>ex)"eX' 
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Step  3.5  requires  Newton-Kaphson  iteration  to  solve  the  six 
nonlinear  equations  for  the  six  unknowns.  In  addition, 
numerical  integration  is  needed  to  calculate  the  generalized 
mass  and  stiffness.  For  the  initial  cycle,  this  is  a  time 
consuming  process  and  will  require  50  to  100  iterations. 
Subsequent  cycles  should  have  sufficiently  small  changes  in 
the  variables  to  require  less  than  ten  iterations. 

STEP  3.6  Determine  amplitude  scaling  factor 


hK2  ♦"(1/2) 

$(x)  =  B  <Kx) ,  scale  eigenvector  accordingly. 


STEP  3.7  Calculate  plastic  zone  dislocation  density 


Ap  =  ~  [1  -  exp(-l/2XEp)]  +  A0exp(-l/2X£p). 

STEP  3.8  Compute  average  amplitude-dependent  damping  ratio 

I?.  =  [A  \  e20„(e )dV  +  A  f  e20„(E)dV]/f  E2dV 

H  1  e  'e  H  P  P  H 

Step  3.8  has  several  numerical  integration  computations  and 
computes  the  exponential  integral  function  inside  an 
integral.  This  is  a  time  consuming  process  and,  along  with 
step  3.5,  will  determine  the  overall  program  runtime. 

STEP  3.9  Compute  frequency-dependent  damping  ratio 

2 

.  2Trot  ET  wt 
4I  "  2  2' 

1  +  W  T 

STEP  3.10  Sum  average  amplitude-  and  frequency-dependent 
damping 


A  =  AH  +  Ar 


STEP  3.11  Output  results  for  (j 
Fatigue  cycle  number,  j 


Fatigue  crack  length,  a 
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Stress  intensity  factor, 

Torsional  stiffness  at  crack, 

Fraction  of  life  expended,  LF 
Plastic  zone  size,  r^ 

Plastic  zone  strain,  t 

P 

Resonant  frequency,  w 

Average  amplitude -dependent  damping  ratio, 
Frequency-dependent  damping  ratio,  Aj 
Total  damping  ratio,  A 
STEP  3.12  Increment  fatigue  cycle  counter 

j++ 

a  ,  =  a  ,  reset  crack  length  variable, 

old  new 

STEP  4.  stop 
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/*  <p>  d:\lc\header\  */ 

j*  ********************************************* ************************ 

NAME 


dmp_crck.exe . calculate  internal  damping  for  a  beam 

undergoing  fatigue  crack  growth 


SYNOPSIS 


dmp_crck  <data  filename>  [data  plot  filename]  execute  program 

dmp_crck  with  input  data  provided  via  a 
data  file 

data  filename  file  containing  input  data  of  beam  and 

material  characteristics 


data  plot  filename 
DESCRIPTION 


file  containing  output  data  for  DISSPIA 
plots  on  PRIME 


The  internal  damping  of  a  beam  undergoing  fatigue  crack  growth  is 
calculated  based  on  an  estimation  of  the  plastic  zone  size  and 
material  dependent  constants  related  to  crack  growth  rate  and 
internal  damping.  The  following  algorithm  is  used  to  calculate 
the  internal  damping  of  a  four-point  beam  with  a  propagating 
crack. 


STEP  1.  Input  material  and  crack  growth  constants 


a_0 
a_disl 
a_f inal 
alpha 
b_crack 

c 

c  1 
c~2 
C_heat 
chi 

delta_ext 

e 

ei 

njnax 

n_ult 

g_check 

growth_rate 


h 

k 

k_knife_edge 

l 

lambda_0 
l_exci ter 

lload 

m_beam 
m_exci ter 
P 

e_yield 

slope_factor 

temperature 


=  initial  crack  length  (in). 

=  dislocation  density  coefficient  (disl/in'2). 

=  final  crack  length  (in). 

=  linear  expansion  coefficient  (in/in/xF). 

=  material  constant  relating  crack  length  to 
torsional  stiffness  (lb- in). 

=  Manson-Coff in  ductility  exponent. 

=  internal  damping  coefficient. 

=  internal  damping  coefficient  (in/in). 

=  specific  heat  per  unit  volume  (psi/xF) 

=  strain  independent  probability  for  the 
remobilization  of  immobile  dislocations. 

=  damping  ratio  due  to  external  source 
=  Young's  modulus  (psi). 

=  Young's  modulus  x  moment  of  inertia  (lb-in'2). 
=  maximum  strain  in  beam  during  danping  data 
collection  (in/in). 

=  static  plastic  fracture  value  at  1/2  cycle 
(in/in). 

=  minimum  crack  growth  before  updating 
eigenvector  results  (in). 

=  crack  growth  rate  constant  based  on 
dislocation  theory  of  crack  growth 
(in'7/lb'A-cycle). 

=  beam  height  (in). 

=  thermal  conductivity  (Ib/sec-xF) 

=  torsional  stiffness  at  knife-edge  support 
(lb-in). 

=  beam  length 

=  dislocatin  density  for  n_p  =  0  (disl/in'2). 

=  location  of  exciter  relative  to  knife-edge 
support  (in). 

=  location  of  load  relative  to  knife-edge 
support  (in). 

=  mass  of  beam  per  unit  length  (slug/ inch). 

=  mass  of  exciter  and  counterweight  (slug). 

=  fatigue  load  applied  to  beam  (lb). 

=  yield  stress  (psi). 

=  constant  relating  beam  slope  and  damping. 

=  ambient  temperature  of  specimen  (xF) 
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STEP  2.  Initialize  variables 


a_old  =  a_0  crack  length 

j  =1  fatigue  cycle  counter 

life  =0.  life  expended  variable 

tau  =  h‘2*c_heat/c‘2*k  relaxation  parameter 

STEP  3.  while  a_last  <=  a_final 

STEP  3.1  Calculate  crack  length  for  current  cycle  based  on 

'a_new 

(da/ki'4)  =  growth_rate. 

.  a_old 

Begin  by  calculating  the  stress  intensity  factor 
using  an  empirical  formula 

f (a/h)  =  1.122  -  1.4*(a/h)  ♦  7.33*(a/h)‘2  - 
13.08*(a/h>*3  ♦  14*(a/h)‘4, 

or 

f(a/h)  =  dsqrt(tan(psi )/psi )* 

(.923  ♦  .199*0  -  sinCpsi ))'4/cos(psi ), 

where  psi  =  c*a/(2*h) 

k_crack  =  b_crack*(h/a) ‘3  torsional  stiffness 

at  crack 

Rotation  of  beam  at  knife-edge  support 

theta  =  (2*ei  +  l_load*k_crack)*l_load*p/ 
<4*ei*(k_knife_edge  ♦  k_crack) 

Nominal  bending  moment  between  load  points 

moment  =  .5*l_load*p  -  k_kni fe_edge*theta 
e_nominal  =  6*moment/h‘2  nominal  bending  stress 
ki  =  e_nominal*dsqrt(c*a)*f(a/h)  stress  intensity 

factor 


STEP  3.1.1 
STEP  3.1.2 
residual 


Assume  a  value  for  a[j] [0] 

Use  Newton-Raphson  Interation  with 
’a_new  ( i  - 1  ] 

=  (da/(ki(a))"4)  -  growth_rate, 

a  old 


slope  =  1/(ki(a_new[i-1]))‘4, 

a_new[i]  =  a_old(i-1]  -  residua l /si ope. 

This  step  requires  nunerical  integration  and  Newton-Raphson 
iteration.  This  case  has  one  equation  and  one  unknown  and  should 
not  require  extensive  computing. 

STEP  3.2  Sum  the  new  value  of  life  expended 

delta_a  =  a_new  -  a_old  incremental  crack  growth 
n_p  =  n_ult*(a_new/a_f inaD'c  plastic  strain 
life  ♦=  2*n_p/n_ult  include  2  only  if  loading  occurs  on 
both  up  and  down  stroke  of  piston 

STEP  3.3  Compute  new  plastic  zone  size 


r_plastic  =  (ki/e_yield)‘2/(c*6). 


where  ki  is  computed  from  a_new. 

STEP  3.4  Calculate  new  torsional  stiffness  at  crack 
k_crack  =  b_crack*(h/a_new)'3 

STEP  3.5  Compute  beam  fundamental  mode  characteristics 
including  eigenvalue,  kappa,  eigenvector,  phi(x),  and 
resonant  frequency,  omega.  The  six  unknowns,  kappa,  four 
from  phi(x)  and  omega,  are  computed  from  the  four  boundary 
conditions: 


w(0)  =  0  deflection  at  knife_edge  is  zero 

ei*w"(0)  =  k_knife_edge*w'(0)  moment  at  knife-edge  is 

equal  to  the  torsional 
stiffness  x  slope 

ei*w"(l/2)  =  k_crack*w'(l/2)  moment  at  crack  is  equal 

to  the  torsional  stiffness 
x  slope 

w,,l(l/2)  =  0  shear  at  crack  is  zero 


Frequency  equation  (ignoring  the  small  amount  of  damping): 


omega‘2*m_general  -  k_general  =  0 


where 


m_general 


■l 

m  beam* phi (x)'2  dx  ♦  phi(l_exciter)'2*m_exciter 

0 


and 


k_general 


‘I 

ei [phi"(x)] ‘2  dx. 

0 


This  step  requires  Newton-Raphson  iteration  to  solve  the  six 
nonlinear  equations  for  the  six  unknowns.  In  addition,  numerical 
integration  is  needed  to  calculate  the  generalized  mass  and 
stiffness.  For  the  initial  cycle,  this  is  a  time  consulting 
process  and  will  require  50  to  100  iterations.  Subsequent  cycles 
should  have  sufficiently  small  changes  in  the  variables  to  require 
less  than  ten  iterations. 


STEP  3.6  Determine  amplitude  scaling  factor 

scale_factor  =  2*n_max/(h*kappa‘2*phi"< 1/2)) 

phi(x)  =  scale_factor*phi(x)  scale  eigenvector  accordingly 

STEP  3.7  Calculate  plastic  zone  dislocation  density 

lambda_p  =  (  (a_disl/chi )*(1-exp( .5*chi*n _p)  ♦ 
dsqrt(lambda_0)*exp< ,5*chi*n_p)  )‘2 

STEP  3.8  Compute  average  amplitude-dependent  damping  ratio 

delta_h  =  (  lambda_p*f  n(x,z)‘2*theta_h(n(x,z))  dV  + 

Jp 

lambda_0*j  n(x,z)‘2*theta_h(n(x,z))  dV  )/ 

J  n(x,z)‘2  dV 

where  the  functions  n(x,z)  and  theta_h(n(x,z))  are  defined  by 
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n(x,z)  =  z*scale_factor*kappa'2*phi"(x) 
theta_h(n(x,z))  =  c_1*(c_2/n(x,z)) ‘2*ei 1 (c_2/n(x, z) ) 

and  ei1()  is  the  exponential  integral  function. 

This  step  has  several  numerical  integration  computations  and 
computes  the  exponential  integral  function  inside  an  integral. 
This  is  a  time  consuming  process  and  along  with  step  3.6  will 
determine  the  overall  program  runtime. 

STEP  3.9  Compute  frequency  dependent  damping  ratio 

delta_i  =  (2*c*alpha‘2*e*temperature)/c_heat* 
omega*tau/(1  ♦  (omega*tau)‘2) 

STEP  3.10  Sum  average  amplitude-  and  frequency  dependent 
damping 

delta  =  detta_i  +  delta_h 

STEP  3.11  Output  results  for  jXlOO  ==  0: 

fatigue  cycle  number,  j 
fatigue  crack  length,  a_new 
stress  intensity  factor,  ki 
torsional  stiffness  at  crack,  k_crack 
fraction  of  life  expended,  life 
plastic  zone  size,  r_plastic 
plastic  zone  strain,  n_p 
resonant  frequency,  omega 

average  ampl i tude- dependent  damping  ratio,  delta_h 
frequency  dependent  damping  ratio,  delta_i 
total  damping  ratio,  delta 

STEP  4.  stop 

FUNCTIONS  USED 


beamjnode . compute  beam  eigenvector,  eigenvalue 

and  frequency 

crack_integral . compute  crack  length  residual 

d  modal . derivative  of  beam  characteristic 

equation  coefficients 

da_dn . change  in  crack  length  as  a  function  of 

fatigue  cycle  number 

data_plot . output  to  file  data  for  DISSPLA  plot 

data_print . output  to  screen  data  printout 

ddphi . second  derivative  of  beam 

characteristic  equation  with  respect 
to  the  beam  coordinate 

dexp . Lattice  SSP  exponential  function 

dsqrt . Lattice  SSP  square  root  function 

ezb_limit . starting  limit  for  second  integration 

of  multiple  integration 

eze_limit . ending  limit  for  second  integration  of 

multiple  integration 

strain2 . strain  squared  function 

n2_theta_h . strain  squared,  amplitude  dependent 

damping  function 

f . stress  intensity  factor  empirical 

formula  for  pure  bending 

get_data . input  material  data  from  file 

input_print . print  input  data  to  screen 

ki . stress  intensity  factor 

mat_assign . MATRIX  structure  assignment  function 

modal . beam  characteristic  equation 

coefficients 

multi_al . multiple  integral  algorithm 

newton_raphson . solves  nonlinear  equation  using  Newton- 
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Raphson  algorithm 

phi . beam  characteristic  equation 

Ph'2 .  . square  of  beam  characteristic  equation 

power . Lattice  SSP  power  function 

pzb_limit . starting  limit  for  second  integration 

of  multiple  integration 

pze_limit . ending  limit  for  second  integration  of 

multiple  integration 

simpson . numerical  integration  using  Simpson's 

rule 
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*********************************************************************  */ 


#include  <stdio.h> 
^include  <math.h> 
#include  <numeric.h> 
#include  <storage.h> 


#define  cube(x)  ((x)*(x)*(x))  /* 

#define  sqr(x)  (<x)*(x)>  /* 

#define  quad(x)  sqr(x)*sqr(x)  /* 

^define  PVERSION  "1.0"  /* 

#def ine  PDATE  "11  AUG  1987"  /* 

/*  program  title  */ 
tfdefine  PTITLE  "Internal  damping  vs 
#def ine  MAX  I  NT  6  /* 

^define  MAX_REAL  36  /* 

unsigned  _STACK  =  4096;  /* 


cube  function  */ 
square  function  */ 
fourth  power  function  V 
program  version  number  */ 
program  date  */ 

crack  propagation  model" 
number  of  integers  input  from  file  */ 
number  of  reals  input  from  file  */ 

stack  size  */ 


«ain(  argc,  argv  ) 
int  argc; 
char  *argv[]; 

double  a_0,  adisl,  a_final,  a_new,  a_old,  b_crack,  c,  c_1,  c_2,  c_heat; 
double  chi,  delta,  delta_crack,  deita_ext,  delta_h,  delta_i,  depth,  e; 
double  ei,  g_check,  g_last,  gen_mass,  growth_rate,  h,  k,  kappa,  k_crack; 
double  k_knife_edge,  l,  lambda_0,  lambda_p,  T_exciter,  life,  l  load; 
double  m_beam,  m_exciter,  moment,  omega,  p,  r_plastic,  sif,  slope_factor; 
double  tau,  temperature,  theta,  toler; 
double  alpha,  njnax,  n_p,  n_ult,  e_nominal,  e_yield; 


/*  double  arrays  */ 

double  crack  19] ,  ez_limit[2],  input_real [MAX_REAl] ,  mode_data [8] ; 
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double  phi_data [5] ,  plot_pak[7],  pr int_pak [12] ,  pz_timit[2],  theta_h[7]; 
double  theta_p[7); 

/*  double  functions  */ 

double  crack_integral(),  da_dn(),  ddphiO,  dexpf),  dphiO,  dsqrtf); 
double  ezb_l imi tf),  eze_limitO,  ff),  kiO,  multi_al(); 
double  newton_raphsonO ,  power(),  pzb_limitf),  pze_limitf ),  simpsonf); 
double  strain?),  strain2(),  n2_theta_h[),  phif),  ph i 2C ) ; 

D_SUPPL0  a_fun,  da_fun,  e_begin,  e_end,  e_fun,  p_begin,  p_end,  p_fun; 
D_SUPPLD  n  begin,  n  end,  n  fun,  h_fun; 

FILE  *f1,  *f2; 

FILE  "fopenf ); 

int  data_rate,  i,  j,  m,  maxim,  n; 

int  i npu t_ i  nt  [HAX_  I N  T  ] ; 

int  data _plot(),  fcloseO,  get  dataO; 

MATRIX  mode; 

MATRIX  beamjnodeO,  d_modal(),  mat_assignf),  modalO; 

M_SUPPLD  m_fun,  dm_fun; 

unsigned  int  is_type_sif,  is_plot_f i le; 

void  data _pr into,  input __pr into; 

printf (  «  Title:  Xs\n",  PTITLE  ); 
printff  "  version:  Xs\n",  PVERSION  ); 

printfC  "  Date:  Xs\n\n",  PDATE  ); 

/*  check  for  data  file  in  command  line  */ 

if(  (  argc  <  2  )  ||  (  argc  >  3  )  )  {  /*  invalid  number  of  data  files  */ 

printf (  11  \n  Usage:  programname  filename  [output  filename]"  ); 
exit(  1  ); 

) 

/*  open  data  file  */ 

iff  (  fl  =  fopenf  argvtl],  »r"  )  )  ==  NULL  )  { 
printff  "  Cannot  open  Xs\n",  argv[1]  ); 
exitf  1  ); 

) 

/*  make  or  overwrite  output  file  */ 

is_plot_file  =  FALSE;  /*  do  not  output  to  file  */ 

if  (  argc  ==  3  )  { 

if  (  (  f2  =  fopenf  argv[21,  "w"  )  )  ==  NULL  >  { 
printff  “  Cannot  open  Xs\n",  argv[2]  ); 
exitf  1  ); 

) 

else  { 

is _plot_file  =  TRUE;  /*  output  to  file  */ 

printff  "  Creating  data  output  file  Xs\n\n",  argv[2]  ); 

> 


> 


/*  STEP  1.  Input  material  and  crsck  growth  constants 

in  addition  to  material  and  crack  growth  constants,  the  following 
program  control  data  must  also  be  input: 


data_  ate  =  number  of  fatigue  cycles  between  data  output 
is_type_sif  =  0  for  constant  load 

1  for  constant  stress  intensity  factor 
j  =  initial  cycle  nutber 

m  =  number  of  subintervals  for  integration 

maxum  =  nurber  of  iteration  cycles 

n  =  nuiber  of  subintervals  for  second  integration 

toler  =  accuracy  tolerance  for  Newton-Raphson  iteration 


*/ 


if  f  get_dataf  fl,  input_int,  input_real  )  )  f 

printff  "  Error  in  data  input  from  Xs\n",  argvtl]  ); 
exitf  1  ); 
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> 


if  (  fclose(  fl  )  )  /*  finished  with  input  file  */ 

( 

printff  "  Error  in  closing  input  data  file  %s\n",  argvll]  ); 
exit(  1  );  /*  could  not  close  data  file  */ 

> 

/*  print  input  data  to  screen  */ 
input_print(  input_int,  input_real  ); 


/*  set  constants  from  input  file  */ 


a_0 

=  input_real [0]  ; 

a_disl 

=  input_real [1]  ; 

a_f inal 

=  input_real [21  ; 

alpha 

=  input_real [33 ; 

b_crack 

=  input_real [4]  ; 

c 

=  input_real [5]  ; 

c  1 

=  input_real [6]  ; 

c~2 

=  input_real [7]  ; 

c_heat 

*  input_real [8]  ; 

chi 

=  input_real [9]  ; 

data_rate 

=  input_int [0]  ; 

delta_ext 

=  input_real [10]  ; 

e 

=  input_real [11]  ; 

ei 

=  input_real [12]  ; 

n_max 

=  input_real [13] ; 

n_ult 

=  input_real [14]  ; 

g_check 

=  input_real [15]  ; 

growth  rate 

=  input_real [16] ; 

h 

=  input_real [17]  ; 

is_type_si f 

=  input_int[1]  ; 

j 

=  input_int[2]; 

k 

=  input_real [18]  ; 

k  kni fe  edge 

=  input_real [19] ; 

r 

=  input_real [20]  ; 

l_exciter 

=  input_real [21] ; 

l_load 

=  input_real [22]  ; 

lambdaO 

=  input_real [23]; 

m 

*  input_int [3] ; 

maxun 

=  input_int [4] ; 

m_beam 

=  input_real [24] ; 

m_exciter 

=  input_real [25] ; 

n 

=  input_int[5]  ; 

P 

=  input_real [26] ; 

e_yield 

=  input_real [27]  ; 

slope_factor 

=  input_real [28]  ; 

temperature 

=  input_real [29] ; 

toler 

=  input_real [30] ; 

/*  STEP  2. 

Initialize  variables 

*/ 

a_old 

=  a  0;  /* 

crack  length  */ 

depth 

=  0.5;  /* 

plastic  zone  depth  factor  V 

g_last 

=  0.;  /* 

eigenvector  update  parameter  */ 

life 

=  0.;  /* 

life  expended  variable  */ 

phi_data[4] 

=  1.;  /* 

assume  an  initial  scale  factor  of 

tau 

-  sqr<  h/PI  )*c  heat/k;  /*  relaxation  constant  */ 

mode 

=  mat_assign(  5,  1, 

NULL  );  /*  set  aside  memory  for 

beam  mode  characteristics  V 

/*  initial  guess  of  mode  variables  */ 
for  (  i  =  0;  i  <  5;  i++  ) 

mode. element [i]  =  input_real [31  ♦  ij; 

a_fun.funct  =  crack_integral;  /*  crack  length  integral  function  pointer  */ 
da_fun.funct  =  da_dn;  /*  da/dN  function  pointer  */ 

a_fun.coeff  =  crack;  /*  the  same  coefficients  are  required  to  */ 

da_fun.coeff  =  crack;  /*  calculate  the  crack  length  and  the  */ 
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crack [0] 

=  b_crack; 

/*  da/dN  function. 

*/ 

cracktl] 

=  ei; 

/*  pass  seven  coefficients  to  a_fun  & 

V 

crack  [21 

=  growth_rate; 

/*  df_fun. 

V 

crack  [3] 

=  h; 

/*  these  coefficents  calculate  e_nominal 

*/ 

crack  [4] 

=  k_kni fe_edge; 

/*  which  is  used  to  determine  the  stress 

*/ 

crackCS] 

=  l_load; 

/*  intensity  factor  on  which  the  functions 

*/ 

crack [6] 

=  p; 

/*  are  based. 

*/ 

crack  [7] 

=  (double)  m; 

/*  number  of  subdivisions  for  integration 

*/ 

m_fun. funct 

=  modal; 

/*  beam  mode  characteristics  equations 

*/ 

dm_fun. funct 

=  djnodal; 

/*  derivative  of  characteristics  equations 

*/ 

mfun.coef f 

=  mode_data; 

/*  the  same  coefficients  are  required  to 

*/ 

dm_fun.coef f 

=  mode_data; 

/*  calculate  the  characterises  equations 

V 

mode_data[0] 

=  k  knife  edge; 

/*  and  their  derivatives,  knife-edge 

*/ 

mode_data[2] 

=  l; 

/*  stiffness,  beam  length 

*/ 

modedata  [3] 

-  m_beam; 

/*  beam  mass  per  unit  length 

*/ 

mode_data[4) 

=  m_exciter; 

/*  mass  of  exciter 

*/ 

mode_data[5] 

=  ei; 

/*  modulus  and  moment  of  inertia 

*/ 

mode_data[6] 

=  l_exciter; 

/*  position  of  exciter  relative  to  knife- 

*/ 

/*  edge 

*/ 

mode_data[7] 

=  (double)  m; 

/*  number  of  subdivisions  for  integration 

*/ 

e_beg in. funct 

=  ezb_limit; 

/*  lower  elastic  limit  function 

*/ 

eend. funct 

=  eze_limit; 

/*  upper  elastic  limit  function 

*/ 

e_fun. funct 

=  n2_theta_h; 

/*  strain  squared  x  amplitude-dependent 

*/ 

e_begin.coef f 

=  ez_limit; 

/*  damping  function 

V 

e_end.coef f 

=  ez_limit; 

e_fun.coef f 

=  theta  h; 

ez_l i mi t [0] 

=  0.; 

/*  lower  elastic  limit  is  0. 

*/ 

ez_l imi  t  [1] 

=  h/2.  ; 

/*  upper  elastic  limit  is  one  half  h 

*/ 

theta  h  [0] 

=  c  1; 

/*  internal  damping  coefficient 

*/ 

theta_h [1] 

=  c_2; 

/*  internal  damping  coefficient 

*/ 

p_beg in. funct 

=  pzblimit; 

/*  lower  plastic  limit  function 

*/ 

p_end. funct 

=  pze_limit; 

/*  upper  plastic  limit  function 

*/ 

p_f un. funct 

=  n2_theta_h; 

/*  strain  squared  x  amplitude-dependent 

*/ 

pbegin.coef f 

=  pzlimit; 

/*  damping  function 

*/ 

pend.coef f 

=  pz_limit; 

pfun.coef f 

=  theta  p; 

theta_p[0] 

=  c  1; 

/*  internal  damping  coefficient 

*/ 

thetapCI] 

=  c_2; 

/*  internal  damping  coefficient 

*/ 

n_begin. funct 

=  ezb  limit; 

/*  use  elastic  functions  for  integrating 

*/ 

nend. funct 

=  eze_limit; 

/*  strain  over  the  volume  of  the  beam 

*/ 

n_fun. funct 

=  strain2; 

/*  strain 

*/ 

nbegin.coef f 

=  ez_limit; 

n_end.coef f 

=  ez_limit; 

n_fun.coef f 

=  phi_data; 

/*  use  coefficients  of  ddphi 

*/ 

h_fun. funct 

=  phi2; 

/*  beam  characteristic  equation 

*/ 

h_fun.coeff 

=  phi  data; 

/*  eigenvector  data 

*/ 

/*  calculate  constant  stress  intensity  factor  */ 

k_crack  =  b_crack*cube(  h/a_old  );  /*  torsional  stiffness  at  crack  */ 

/*  Rotation  of  beam  at  knife-edge  support  */ 
theta  =  (  2.*ei  +  l_load*k_crack  )*l_load*p/ 

(  4.*ei*(  k_knife_edge  +  k_cracx  )  ); 

/*  Nominal  bending  moment  between  load  points  */ 

moment  =  ,5*l_(oad*p  -  k_kni fe_edge*theta; 

e_nominal  =  6.*moment/sqr(h);  /*  nominal  bending  stress  */ 

sif  =  ki(  a_old,  h,  e_nominal  );  /*  constant  stress  intensity 

factor  */ 

r_plastic  =  sqr(  sif/e_yield  )/(Pl*6.);  /*  constant  plastic  zone 

radius  */ 

/*  STEP  3.  while  a_old  <=  a_final  */ 
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while  (  a_old  <=  a_final  )  { 

/*  STEP  3.1  Calculate  crack  length  for  current  cycle  */ 

/*  STEP  3.1.1  Assume  a  value  for  a_new  */ 
anew  =  a_old  +  growth_rate*quad(  sif  ); 

if  (  !(  is_type_sif  )  ) 

< 

/*  constant  fatigue  load  */ 

/*  STEP  3.1.2  Use  Newton-Raphson  Interation  */ 

crack  [8]  =  a_old;  /*  lower  limit  of  integral  */ 

a_new  =  r,ewton_raphson<  a_new,  maxum,  toler,  a_fun,  da_fun  ); 

> 

/*  STEP  3.2  Sum  the  new  value  of  life  expended  */ 

n_p  =  n_ult*power(  (  a_new  -  a_old  )/a_final,  c  );  /*  plastic 

strain  */ 

life  +=  power<  n_p/n_ult,  l./c  );  /*  include  2  only  if  loading 

occurs  on  both  up  and  dow- 
stroke  of  piston  */ 

/*  STEP  3.3  Calculate  new  torsional  stiffness  at  crack  */ 
k_crack  =  b_crack*cube(  h/a_new  ); 

mode_data[1]  =  k_crack;  /*  pass  on  crack  stiffness  */ 

if  (  is_type_sif  ) 

{ 

/*  calculate  new  load  */ 

e_nominal  =  sif/C  dsqrtf  Pl*a_neu  )*f(  a_new/h  )  >;  /*  nominal 

bending 
stress  */ 

moment  =  sqr(  h  )*e_nominal/6. ;  /*  nominal  bending  moment  */ 

p  =  4.*ei*(  kcrack  +  k_knife_edge  )*moment/ 

(  (  2.*ei  •  l  load*k_kni fe_edge  )*l_load*k_crack  ); 

> 

else 

C 

/*  for  constant  load,  calculate  new  plastic  zone  radius 
Rotation  of  beam  at  knife-edge  support  */ 
theta  =  <  2.*ei  ♦  l_load*k_crack  )*l_load*p/ 

(  4.*ei*(  k_knife_edge  ♦  k_crack  )  ); 

/*  Nominal  bending  moment  between  load  points  */ 

moment  =  .5*l_load*p  -  kkni f e_edge*theta; 

e_nominal  =  6.*moment/sqr(  h  );  /*  nominal  bending  stress  */ 

/*  STEP  3.4  Compute  new  plastic  zone  size  */ 

rplastic  =  sqr(  (  sif  =  ki(  a_new,  h,  e_nominal  )  )  /e_yield  )/ 

( P I *6 . ) ; 

> 

/*  STEP  3.5  Compute  beam  fundamental  mode  characteristics 
including  eigenvalue,  kappa,  eigenvector,  phi(x),  and  resonant 
frequency,  omega.  Use  the  nonlinear  system  of  equations  function 
to  solve  this  5  equation  system.  Assign  the  vector  x  as  follows 

element [0]  -  eigenvalue,  kappa 
element [1]  -  resonant  frequency,  omega 
elemental  -  B  coefficient 
element [3]  -  C  coefficient 
elemental  -  0  coefficient 

where  phi(x)  =  cos(  kappa*x  )  +  B*sin(  kappa*x  )  + 

C*cosh(  kappa*x  )  +  D*sinh(  kappa*x  ), 
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and  ddphi(x)  =  sqr<  kappa  )*(  -cos(  kappa*x  )  -  B*sin(  kappa*x  )  + 
C*cosh(  kappa*x  )  +  D*sinh<  kappa'x  )  )  */ 

if  (  abs(  a_neu  -  g_last  )  >  g_check  )  /*  only  update  for  large 

crack  growth  */ 

{ 

node  =  beam_mode(  mode,  maxum,  toler,  m_fun,  dm_fun  ); 
g_last  =  a_new;  /*  reset  crack  length  of  last  update  */ 

kappa  =  mode. element [0] ;  /*  eigenvalue  V 

omega  =  mode. element [1] ;  /*  resonant  frequency  */ 

/*  STEP  3.6  Determine  anplitude  scaling  factor  */ 

phi_data[0]  =  kappa; 

phi_data[1]  =  mode. element  12] ; 

phi_datat2]  =  mode. element [3] ; 

phi_dataC3]  =  mode. element  14] ; 

phi_data[4]  =  1.; 

phi_dataC4]  *  n_max/strain(  1,  ,5*h,  phi_data  ); 

/*  STEP  3.7  Calculate  plastic  zone  dislocation  density  */ 
lambda _p  =  (  a_disl/chi  )*(  1  -  dexp(  -.5*chi*n_p  )  )  + 
dsqrtl  lambda_0  )*dexp(  -.5*chi*n_p  ); 
lambdap  =  sqr(  lambda_p  ); 

/*  STEP  3.8  Compute  average  ampl i tude-deperdent  damping  ratio 

delta_h  =  (  lambda_o*t  n(x,z)'2*theta_h<n(x,z))  dv  + 
up 

lambda_0*t  n(x, z) '2*theta_h(n(x, z))  dV  )/ 
ue 

t  n(x,z)‘2  dV 
u 

where  the  functions  n(x,z)  and  t  uta_h(nlx,z>)  are  defined  by 

n(x,z)  =  z*scal e_iactor*kappa'2*h"(x) 
theta_h(n(x,z))  =  c_1*(c_2/n(x,z)) '2*ei 1 (c_2/n(x, z) ) 

and  e i  1  < )  is  the  expot.tiu  . -i  >nteg-  l  ^unction.  */ 

/*  check  size  of  plastic  zone  relative  to  remaining  section  depth  */ 
if  (  <  pz_timit[0]  =  0.5*(h  •  a_new)  -  r_plastic  )  <  0.0  ) 

/*  rplastic  is  too  large  for  remaining  section  depth  */ 

{ 

r_plastic  =  0.5*(h  -  a_new);  /*  r_plastic  is  half  of 

remaining  section  depth  */ 


pz  limit  [0] 

=  0.0; 

/* 

treat  plastic  zone  as  full 

depth  of  section 

*/ 

depth 

\ 

=  1.0; 

j 

pz_l imi t [1] 

=  0.5*<h  -  a_new); 

r 

upper  plastic  limit  is  one 

half  h  less  crack  length 

*/ 

theta_h [2] 

= 

kappa; 

/* 

eigenvalue 

V 

theta_h [3] 

= 

mode. element  C2] ; 

!* 

phi<x)  coefficient  B 

*/ 

thela_h  14] 

= 

mode. element  13] ; 

r 

phi(x)  coefficient  C 

*/ 

theta_h [5] 

= 

mode. element [4] ; 

r 

phi(x)  coefficient  D 

*/ 

theta_h (6] 

= 

phi_data[4] ; 

/* 

amplitude  scale  factor 

*/ 

theta_p[2] 

= 

kappa; 

/* 

eigenvalue 

V 

theta_p[3] 

= 

mode. element [2] ; 

/* 

phi(x)  coefficient  B 

*/ 

theta_pl4] 

= 

mode. element [3] ; 

/* 

phi(x)  coefficient  C 

*/ 

theta_p[5] 

= 

mode. element [4] ; 

/* 

phi(x)  coefficient  D 

*/ 

theta_p[6) 

= 

.25*sqr(  h/pz  limitll] 

)*phi  data[4]; 

/* 

amplitude  scale  factor 
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adjusted  for  reduced  cross- 
section  in  plastic  zone  V 

del ta_h  =  (  multi_al(  0.,  I,  m,  n,  e_begin,  e_end,  e_fun  )  ♦ 
depth*(  (  lambda_p  -  lambda_0  )/lambda_0  )* 
multi_al(  t  -  r_plastic,  l,  m,  n,  p_begin,  p_end,  p_run  )  )/ 
multi_al(  0.,  I,  m,  n,  n_begin,  n_end,  n_fun  ); 

/*  STEP  3.9  Compute  frequency  dependent  damping  ratio  */ 

delta_i  =  (2.*Pl*sqr(  alpha  )*e*temperature)/c_heat* 
omega*tau/(  1  +  sqr(  omega*tau  )  ); 

/*  calculate  damping  due  to  crack  surfaces  based  on  slope  of  beam  at 
crack  */ 

delta_crack  =  slope_factor*abs(  dphi(  l,  phi_data  )  ); 

/*  STEP  3.10  Sun  average  amplitude-  and  frequency  dependent 
dauiping.  Also,  sum  external  and  crack  related  damping  V 

delta  =  delta_crack  +  delta_ext  +  delta_i  ♦  delta_h; 

>  /*  end  of  update  block  */ 

/*  STEP  3.11  Output  results  for  j  X  data_rate  ==  0: 

fatigue  cycle  number,  j 
fatigue  crack  length,  a_new 
stress  intensity  factor,  ki 
torsional  stiffness  at  crack,  k_crack 
fraction  of  life  expended,  life 
plastic  zone  size,  r_plastie 
plastic  zone  strain,  n_p 
resonant  frequency,  omega 

average  amplitude -dependent  damping  ratio,  delta_h 
frequency  dependent  damping  ratio,  delta_i 
total  damping  ratio,  delta  */ 

if  (  (  j  X  data_rate  )  ==  0  )  < 

if  (  is_plot_file  )  < 


plot_pak [0]  =  p;  /*  fatigie  load  */ 
plot_pak[1]  =  e_nominal/e;  /*  fatigue  nominal  strain  */ 
plot_pak[2]  =  a_neu;  /*  fatigue  crack  length  */ 
plot_pak[3]  =  sif;  /*  stress  intensity  factor  */ 
plot_pak[4j  =  delta;  /*  total  internal  damping  */ 


plot_pak[5]  =  omega/(2.*Pl );  /*  resonant  frequency  (hz)  */ 

/*  compute  ratio  of  max  beam  strain  to  input  forcing  function, 
generalized  mass  */ 

genjnass  =  m_beam* simps on (  0.,  I,  m,  h_fun  )  + 

m_exciter*phi2(  l_exciter,  phi_data  ); 
plot_pak[6]  =  n_max*abs(  phi<  l_exciter,  phi_data  )  )/ 

(  2.*delta*sqr(  omega*gen__mass  )  ); 
if  (  data_plot(  f2,  j,  plot_pak  )  !=  0  5  /*  output  results  in 

plot  format  */ 
t 

printft  "Fite  write  error  in  main\n"  ); 
exit(  1  ); 

> 

>  else  C 


print_pak[0]  =  a_new;  /*  fatigue  crack  length  */ 
print_pak[1J  =  sif;  /*  stress  intensity  factor  */ 
print_pak[2]  =  k_crack;  /*  torsional  stiffness  at  crack  */ 
print_pak[3]  =  life;  /*  fraction  of  life  expended  */ 
print_pak[4J  =  r_plastic;  /*  plastic  zone  size  */ 
print_pak[5]  =  n_p;  /*  plastic  zone  strain  */ 
print_pak[6]  =  omega;  /*  resonant  frequency  */ 
print_pakl7J  =  delta_h;  /*  amplitude-dependent  damping  */ 
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print_pak[8]  =  delta_i;  /*  frequency  dependent  damping  */ 


print _pak[9]  =  delta_crack;  /*  crack  dependent  damping  */ 
print _pak[10]  =  delta;  /*  total  damping  ratio  */ 
print_pak [1 1]  =  p;  /*  fatigue  load  */ 
data_print(  j,  print  pak  );  /*  put  results  on  screen  */ 


> 

> 

a_old  =  a_new;  /*  increment  crack  length  */ 

♦+j;  /*  increment  cycle  counter  */ 

>  /*  end  of  while  loop  */ 

i f  (  is _plot_f i le  ) 

if  (  fclosel  f?  )  )  /*  finished  with  output  file  */ 

{ 

printfl  11  Error  in  closing  output  data  file  Xs\n'*(  argv[2]  ); 
exitC  1  );  /*  could  not  close  output  file  */ 

> 

/*  STEP  4.  Stop  */ 

printfC  "  Normal  exit\n"  ); 
exi t(  0  ); 
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/*  <p>  d:\lc\header\  */ 

/*  ********************************************************************* 


This  set  of  routines  provides  specific  support  for  the  program 
dmp_crck.  These  functions  are  for  use  with  dmp_crck  only  and  may 
require  modification  to  work  with  other  programs. 


beam_mode . compute  beam  eigenvector,  eigenvalue 

and  frequency 

crack_integral . compute  crack  length  residual 

d_modal . derivative  of  beam  characteristic 

equation  coefficients 

dadn . change  in  crack  length  as  a  function  of 

fatigue  cycle  number 

data_plot . output  to  file  data  for  DISSPLA  plot 

data  print . output  to  screen  data  printout 

dddphi . third  derivative  of  beam  characteristic 

equation  with  respect  to  the  beam 
coordinate 

ddphi . second  derivative  of  beam 

characteristic  equation  with  respect 
to  the  beam  coordinate 

ddphi 2 . square  of  second  derivative  of  beam 

characteristic  equation  with  respect 
to  the  beam  coordinate 

dphi . first  derivative  of  beam  characteristic 

equation  with  respect  to  the  beam 
coordinate 

ezb_limit . starting  limit  for  second  integration 

of  multiple  integration 

eze_limit . ending  limit  for  second  integration  of 

multiple  integration 

f . empirical  formula  for  stress  intensity 

factor  in  pure  bending 

get_data . input  material  data  from  file 

input  print . print  input  data  on  screen 

ki . stress  intensity  factor  for  pure 

bending 

modal . beam  characteristic  equation 

coefficients 

strain2 . strain  squared  function 

n2_theta_h . strain  squared,  amplitude  dependent 

damping  function 

phi . beam  characteristic  equation 

phi2 . square  of  beam  characteristic  equation 

pzb_limit . starting  limit  for  second  integration 

of  multiple  integration 

pze_limit . ending  limit  for  second  integration  of 

multiple  integration 

strain . beam  strain  function 


ft********************************************************************  *  j 


#include  <stdio.h> 
#include  <math.h> 
^include  <numeric.h> 
^include  <storage.h> 


#define  cube(x)  (<x)*(x)*(x))  /* 
tfdefine  sqr(x)  ((x)*(x))  /* 
/(define  quad(x)  sqr(x)*sqr(x)  /* 
#de*''ne  MAX_INT  6  /* 
#def inc  .!AX_PFAl  36  /* 


cube  function  */ 

square  function  */ 

fourth  power  function  */ 

number  of  integers  input  from  file  */ 

number  of  reals  input  from  file  */ 
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MATRIX  beamjaodeC  mode,  maxum,  toler,  mode_fun,  dmode_fun) 

I*  ********************************************************************* 

NAME 

beam_mode . Compute  beam  fundametal  characteristics 

including  eigenvalue,  eigenvector,  and 
resonant  frequency 


SYNOPSIS 


include  <numeric.h> 


mode  =  beam_mode(  mode,  maxixn,  toler,  mode_fun,  dmode_fun);  beam 

characteristics  computed  by  Newton- 
Raphson  iteration 


MATRIX  mode; 

int  maxum; 

double  toler; 
M_SUPPLD  modefun, 

dmode_fun; 


vector  containing  initial  guess  of 
characteristics  as  input,  and  final 
solution  values  as  output 
maximum  number  of  Newton- Raphson 
iterations 

accuracy  tolerance  on  final  solution 
user  supplied  beam  characteristic 
equations  and  coefficients 
user  supplied  derivative  of  beam 
characteristic  equations  and 
coefficients 


DESCRIPTION 

A  vibrating  beam  has  four  boundary  conditions  and  an  equation  of 
motion.  Applying  the  Newton-Raphson  iteration  algorithm  to  these 
five  equations  results  in  a  solution  for  the  eigenvalue, 
eigenvector  and  resonant  frequency. 

FUNCTIONS  USED 


REFERENCES 

Burden,  R.L.  and  Faires,  J.D.  "Numerical  Analysis,”  third 
edition,  Prindle,  Weber  and  Schmidt,  Boston,  1984,  pp.  153- 
195. 

Harris,  C.M.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  W.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  1987 


V 


double  toler; 
int  maxum; 

MATRIX  mode; 

M  SUPPLD  dmodefun,  mode  fun; 

“  < 

MATRIX  newtonO; 

/*  begin  newton- raphson  iterative  solution  */ 

mode  =  newton(  mode,  maxim,  toler,  mode_fun,  dnode_fun  ); 

/*  check  for  convergence  */ 
if  (  mat_err  ==  TRUE  ) 

printff  "newton  solution  failed  to  converge  \n"  >; 
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/*  return  vector  solution  */ 
return  (  mode  ); 

} 


dociile  crack  Integra l (  a,  co  ) 

I*  ........................... 

NAME 


crack_integral . computes  crack  length  residual  integral 

using  Simpson's  Rule 

SYNOPSIS 

include  <nuneric.h> 


delta_a  =  crack_integral(  a,  co  );  Compute  crack  length  residual 

from  integral  of  da/dN  relation 


double  a, 

*co, 

delta_a; 


crack  length  guess 

array  of  coefficients  for  use  in  function 
residual  due  to  error  in  estimate  of  a 


DESCRIPTION 


This  function  uses  Simpson's  rule  to  calculate  an  increase  in  the 
crack  length  integral.  The  function  then  computes  the  error  in  the 
crack  length  integral  increase  by  subtracting  the  known  solution 
from  the  numerical  integration  value.  The  resulting  error  is 
returned  to  the  calling  function. 


ta  new 

residual 

=  3  ~  (da/(k 

ua_old 

co  [0] 

=  b_crack 

co[1] 

=  ei 

co  [2] 

=  growth  rate 

co[3] 

=  h 

co  [4] 

=  k_knife_edge 

co  [5] 

=  lload 

cot6] 

*  P 

co  [7] 

=  (double)  m 

co[8] 

=  a  old 

growth_rate. 


FUNCTIONS  USED 


da_dn . change  in  crack  length  at  a  function  of 

fatigue  cycle  number 

simpson . numerical  integration  using  Simpson's 

rule 


REFERENCES 

Frost,  N.E.,  Harsh,  K.J.  and  Pook,  L.P.  "Hetal  Fatigue," 
Clarendon  Press,  Oxford,  1974,  pp  234-235. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  1987 

***********************************.*******************.****.********  * j 

double  a,  *co; 

< 

double  da_dn(),  simpsonO; 

D_SUPPID  da_fun; 
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/*  da/dN  function  pointer  */ 
/*  da/dN  array  pointer  */ 


dafun.f unc t  =  da_dn ; 
da_fun.coeff  =  co; 


/*  return  integral  -  crack  growth  constant  */ 

return  (  simpson(  co[8],  a,  (  int  )  coC7],  da_fun  )  -  1.  ); 


MATRIX  d_»odal(  x,  co,  df  ) 

I*  ********************************************************************* 

NAME 

djnodal . confutes  the  matrix  of  coefficient 

derivatives  for  Newton_Raphson 
iteration 

SYNOPSIS 

include  <numeric.h> 

df  =  d_modal<  co,  df  );  compute  derivative  matrix  of  nonlinear 

equation  coefficients 

double  *co;  pointer  to  array  of  coefficient 

parameters 

MATRIX  x,  vector  of  solution  values 

dx;  matrix  containing  derivative  terms  of 

nonlinear  equation  coefficients 

DESCRIPTION 

User  supplied  routine  to  calculate  the  partial  derivatives  of  the 
nonlinear  system  of  equations  as  a  function  of  x  with  modifiable 
coefficients 

co[0]  =  k_knife_edge; 
co[1]  =  k_crack; 
co[2J  =  lj 
co[3]  =  m  beam; 
coC4]  =  mexciter; 
co  [5]  =  ei; 
co[6]  =  l_exciter; 

co [7]  =  (double)  m; 

FUNCT IONS  USED 

coef f_red. 
dcoef f . . . 
d_determ. . 

dcos . 

dcosh . 

dsin . 

ds i nh . 

mat_assign 
matdelete 

REFERENCES 

Burden,  R.L.  and  Faires,  J.D.  "Numerical  Analysis,"  third 
edition,  Prindle,  Weber  and  Schmidt,  Boston,  1984,  pp.  153- 
195. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
OATE:  18  AUG  87 


reduction  of  coefficients 

partial  derivative  of  coefficient  reduction 

partial  derivative  of  determinate 

Lattice  SSP  cosine  function 

Lattice  SSP  hyperbolic  cosine  function 

Lattice  SSP  sine  function 

Lattice  SSP  hyperbolic  sine  function 

assign  matrix  structure 

delete  matrix  storage  space 


knife-edge  stiffness 
crack  stiffness 
beam  length 

beam  mass  per  unit  length 

mass  of  exciter 

modulus  and  moment  of  inertia 

position  of  exciter  relative  to  knife- 

edge 

number  of  subdivisions  for  integration 
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it********************************************************************  * j 

MATRIX  df,  x; 
double  *co; 

{ 

double  ch_ei,  c_ei,  eigen,  gen_mass,  gen_stiff; 
double  sh_ei,  s_ei,  stiff_c,  stiff_ke; 
double  phi_data[5]; 

double  d_determO,  ddphi2(),  p2k<),  p2kb(),  p2kcO,  p2kd(),  p2m(); 
double  p2mb(),  p2mc(),  p2md(),  ph i 2( > ; 
double  dcosO,  dcoshO,  dsin<),  dsinhO; 
i  nt  i ; 

D_SUPPID  ddh_fun[5] ,  h_fun!5]; 

MATRIX  a,  a_p,  coeff,  coeff_p; 

MATRIX  coeff_red(),  d_coeff(),  mat_assign( ),  mat_delete( ); 

/*  user  supplied  equations  for  beam  analysis  */ 

eigen  =  x. element [0] *co[2] ;  /*  eigenvalue  *  length  */ 

ch_ei  =  dcosh(  eigen  ); 

c_ei  =  dcos(  eigen  ); 

shei  =  dsinh!  eigen  ); 

s_ei  =  dsin(  eigen  ); 

stiff_c  =  co[1]/(  x. element [0] *co!5]  );  /*  characteristic  stiffness 

at  crack  */ 

stiff_ke  =  co[0]/(  x. element  [0]*co[5]  );  /*  characteristic  stiffness 

at  boundary  */ 

/*  characteristic  equation  coefficients  */ 
phi_data[0]  =  x. element [0] ; 

phi|  datatl]  =  x. element [2] ; 

phi_data[2]  =  x. element [3] ; 

phi_data!3]  =  x. elements]  ; 

phi_  data  [A]  =1.;  /*  use  scale  factor  of  1  */ 

/*  create  a  matrix  and  partial  derivative  matrices  */ 
a  =  mat_assign(  4,  4,  NULL  ); 

a_p  =  matassign!  4,  4,  NULL  ); 

/*  fill  a  matrix  */ 
a. element C  0]  =  1  ; 

a. element!  1]  =  1.; 

a.elementC  2]  =  -(  c_ei  ♦  stiff_c*s_ei  ); 

a. element!  3]  =  s_ei; 

a. element!  4]  =  0.; 

a. element!  5)  =  stiff_ke; 

a. element!  6]  =  -s_ei  +  st i f f_c*c_ei ; 

a. element  I  7]  =  -c_ei; 

a. element!  8]  =  1.; 

a. element  19]  =  - 1 .  ; 

a. element  HO]  =  ch_ei  +  st i f f_c*sh_ei ; 

a. element [1 1]  =  sh_ei; 

a. element  112]  =  0.; 

a. element !13]  =  stiff_ke; 

a. element  114]  =  sh_ei  ♦  sti f f_c*ch_ei ; 

a. element [15]  =  ch_ei; 

/*  partials  with  respect  to  eigenvalue  */ 
aja. element!  0]  =  0.; 

a_p. element!  1]  =  0.; 

a_p. element!  2]  =  co[2]*(  (  1.  ♦  st i f f_c/eigen  )*s_ei  -  stiff_c*c_ei  ); 

a _p. element!  3)  =  co!2]*c_ei; 

a_p. element!  4]  =  0.; 

a_p. element!  5]  =  -co[2]*stiff_ke/eigen; 

a_p. element!  61  =  -co!2]*(  (  1.  ♦  st iff_c/eigen  )*c_ei  +  stiff_c*s_ei  ); 
a_p. element!  7]  =  co!2]*s_ei; 
a_p. element!  8]  =0.; 
a_p. element!  9]  =  0.; 
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a _p. element [10]  =  co[2J*(  (  1.  -  stiff_c/eigen  )*sh_ei  ♦  stiff_e*ch_ei  ); 

a_p. element  Ml]  =  co[2]*ch_ei; 

a_p. element [12]  =  0.; 

a_p. element [13]  =  -co[2]*stiff_ke/eigen; 

a_p. element [14]  =  co[2]*(  (  1.  -  stiff_c/eigen  )*ch_ei  ♦  stiff_c*sh_ei  ); 
a  _p. element [15]  =  co[2]*sh_ei; 

/*  0  partial  differential  matrix  */ 
for  (  i  =  0;  i  <=  24;  ++i  ) 
df .element [i]  =  0.; 

/*  partial  derivatives  of  first  equation  with  respect  to  eigenvalue. 

all  other  partials  of  first  equation  are  0  */ 
df .element [0]  =  d_determ<  a,  a_p  ); 

/*  get  coefficients  and  partials  */ 
coeff  =  mat_assign(  4,  1,  NULL  ); 

coeff _p  =  mat_assign(  4,  1,  NULL  ); 

/*  fill  coefficients  in  terms  of  coeff [0]  */ 
coeff  =  coeff_red(  a,  coeff  ); 
coeff _p  =  d_coe?f(  a,  a_p,  coeff_p  ); 

/*  setup  array  of  functions  for  partial  derivatives  of  generalized  mass 
and  stiffness  */ 

h_fun[0] .funct  =  p2m;  /*  mass  with  respect  to  eigenvalue  */ 

ddh_fun[0] . funct  =  p2k;  /*  stiffness  with  respect  to  eigenvalue  */ 

h_fun[1] .funct  =  phi2;  /*  mass  with  respect  to  frequency  */ 

h_fun[2] .funct  =  p2mb;  /*  mass  with  respect  to  coefficient  B  */ 

ddh_fun[2] .funct  =  p2kb;  /*  stiffness  with  respect  to  eigenvalue  B  */ 

h_fun[3] -funct  =  p2mc;  /*  mass  with  respect  to  coefficient  C  */ 

ddh_fun[3] .funct  =  p2kc;  /*  stiffness  with  respect  to  coefficient  C  */ 

h_fun[4] .funct  =  p2md;  /*  mass  with  respect  to  coefficient  D  */ 

ddh_fun[4] .funct  =  p2kd;  /*  stiffness  with  respect  to  coefficient  D  */ 

/*  pass  coefficients  to  integration  functions  */ 
for  (  i  =  0;  i  <  5;  i++  ) 

( 

h_fun[i] .coeff  =  phi_data; 
ddh_fun[i) .coeff  =  phi_data; 

/*  pa^tiat  of  beam  mass  */ 

genjnass  =  co[3]*simpson(  0.,  co[2],  (int)  co[7] ,  h_fun[i]  )  + 

co[4]*(*h_fun[i]  .functK  co[6],  phi_data  ); 

if  (  i  ==  1  ) 

/*  partial  derivative  of  equation  2  with  respect  to  frequency  */ 
df .element [6]  =  2. *x. element [1] *gen_mass; 

else 

{ 

/*  partial  of  beam  stiffness  */ 

gen_stiff  =  co[5]*simpson(  0.,  co[2],  (int)  co[7j,  ddh_fun[i]  ); 
/*  partial  derivatives  */ 

df.elementM  +  5*i]  =  sqr(  x.element[1]  )*gen_mass  -  gen  stiff; 

> 

> 

/*  third  -  fifth  equations,  derivative  with  respect  to  eigenvalue  */ 
for  (  i  -  1;  i  <=  3;  ++i  )  ( 

df .element [i+1]  =  coeff  p. element [i ] ; 

> 

/*  third  -  fifth  equations,  mode  shape  coefficients  */ 
df. element [12]  =  -1.; 
df .element [18]  =  -1.; 
df .element [24]  =  -1.; 
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/*  release  storage  */ 
a  =  mat_delete( 

a_p  =  mat_delete( 

coeff  =  mat_delete( 

coeff_p  =  mat_delete( 


a  ); 
a_p  ); 
coeff  ); 
coeff _p  ); 


return  {  df  ); 

> 


double  da_di(  a,  co  ) 
/*  ****************** 
NAME 


da  dn . 

...crack 

growth  per  fatigue  cycle 

SYNOPSIS 

rate  =  da_dn(  a,  co  ); 

crack 

growth  rate  per  fatigue  cycle 

double  a. 

crack 

length  estimate 

*co, 

array 

of  material  coefficients 

rate; 

crack 

growth  rate 

DESCRIPTION 

Compute  crack  growth  rate  based  on  dislocation  theory.  The 
relation  is 

da/dN  =  alpha*KI‘4 

where  KI  is  the  stress  intensity  factor  and  alpha  is  a  material 
dependent  constant.  The  function  ki  is  used  in  the  computation 
of  the  stress  intensity  factor,  ki  requires  the  following 
calculations 


Calculate  new  torsional  stiffness  at  crack 
k_crack  =  b_crack*cube(  (  h/a_new  )  ); 


Rotation  of  beam  at  knife-edge  support 
theta  =  (  2.*ei  ♦  l_load*k_crack  )*l_load*p/ 
(  4.*ei*(  k_knife_edge  +  k_crack  )  ); 


Nominal  bending  moment  between  load  points 
moment  =  ,5*l_load*p  -  k_kni fe_edge* theta; 
e_nominal  -  6.*moment/sqr(h);  nominal  bending  stress 


coCO] 
coM] 
coC2] 
co  13] 
colA] 
cot5] 
col6J 
co  17] 
co  [8] 


=  b_crack 
=  ei 

=  growth  rate 
=  h 

=  k_knife_edge 
=  l_load 
=  P_ 

=  (double)  m 
=  a  old 


Express  the  relation  in  the  form 
da/alpha*KI '4  =  dN 


FUNCTIONS  USED 


ki . compute  stress  intensity  factor 

REFERENCES 
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Frost,  N.E.,  Marsh,  K.J.  and  Pook,  L.P.  "Metal  Fatigue," 

Clarendon  Press,  Oxford,  1974,  pp  234-235. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  1987 

*********************************************************************  * / 
double  a,  *co; 

double  k_crack,  moment,  e_nominal,  theta; 
double  ki(>; 

/*  Calculate  new  torsional  stiffness  at  crack  */ 
k_crack  =  co[0]*cube(  co[3]/a  ); 

/*  Rotation  of  beam  at  knife-edge  support  */ 
theta  =  (  2.*co[1]  +  co[5]*k_crack  )*co[5)*co[63/ 

(  4.*coC1]*(  co[4]  ♦  k_crack  )  ); 

/*  Nominal  bending  moment  between  load  points  */ 
moment  =  .5*co[5]*co[6]  -  co[4)*theta; 

e_nominal  =  6.*moment/sqr(  co[3]  );  /*  nominal  bending  stress  */ 

/*  return  growth  rate  */ 

return  (  1./(  co[2]*quad(  ki(  a,  co[3],  e_nominal  )  )  )  ); 


int  data_plot(  fp,  n,  data  1 

!+  *************************  ******************************************* 
NAME 


data_plot.  . Output  to  a  file  data  for  use  with 

DISSPLA  programs  on  the  PRIME 

SYNOPSIS 

data_plot(  file,  n,  data  )  Output  data  to  file  in  specific  format 

FILE  *fp;  output  file 

int  n;  fatigue  cycle  number 

double  ‘data;  pointer  to  array  of  output  data 

DESCRIPTION 

This  function  outputs  data  in  the  following  format 

n  p  strain  a_new  ki 

5£d5  %f6  %e11.4  xxxxxxxxxxx  Sell. 4  xxxxxxxxxxx  Sell. 4  xxxxxxxxxxx 


delta  frequency  amplification 
Sell. 4  Sf7. 1  Sell. 4 

where 


ampl i f icat ion 

a_new 

delta 

frequency 

ki 

n 

P 

strain 


=  ratio  of  output  strain  in  beam  to  input  force 
of  exciter  (in/in/lb) 

=  crack  length  (in) 

=  internal  damping  ratio 
=  resonant  frequency  of  beam  (cps) 

=  stress  intensity  factor  (psi-in'1/2) 

=  fatigue  cycle  number 
=  fatigue  load 

=  nominal  strain  in  beam  under  fatigue  load 


and  small  x  represents  a  blank  space. 
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FUNCTIONS  USED 


fprintf . write  to  file 

REFERENCES 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  1987 

*********************************************************************  * j 


double  *data; 

FILE  *fp; 
int  n; 

( 

char  *record  1  =  •'%5d546.0f%11.4E  %11.4E  %11.4E\n"; 

char  *record~2  =  "%1 1 ,4E%7. 1 f%1 1 ,4E\n"; 
int  err; 

/*  initialize  error  code  */ 
err  =  0; 

/*  display  fatigue  cycle  number  and  crack  length  on  screen  */ 
printff  "fatigue  cycle  number:  %5d\n",  n  ); 
printff  "crack  length:  %11.4e  (in)\n“,  data[2]  ); 

printf(  "total  damping:  %11.4e\n",  data[4]  ); 

/*  write  first  record  to  file  */ 

if  (  fprintff  fp,  recordl ,  n,  data[0],  datall],  data[2],  data[3]  )  !=  67 

t 

printff  "File  write  error  in  data_plot\n"  );  /*  error  in  output  */ 

err  =1;  /*  set  error  code  */ 

> 

else  /*  write  second  record  to  file  */ 

if  (  fprintff  fp,  record_2,  data[4],  data[5],  data[6]  )  ! =  30  > 

printft  "File  write  error  in  data_plot\n"  );  /*  error  in  output  */ 
err  =  1;  /*  set  error  code  */ 

> 

return(  err  );  /*  return  error  code  */ 

> 


void  data_print(  n,  data  ) 

I  *  ********************************************************************* 

NAME 

data  print . Output  dmp_crack  results  to  screen 

SYNOPSIS 

data_print(  n,  data  );  Output  crack  growth  and  internal 

damping  results  to  console 

double  *data;  array  of  results 

n;  fatigue  cycle  number 

DESCRIPTION 

Output  table  of  crack  growth  and  internal  damping  results  to 
console.  Table  of  results 

fatigue  cycle  number,  j 
fatigue  crack  length,  a_new 
stress  intensity  factor,  ki 
torsional  stiffness  at  crack,  k_crack 
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fraction  of  life  expended,  life 
plastic  zone  size,  r_plastic 
plastic  zone  strain,  n_p 
resonant  frequency,  omega 

average  amplitude-dependent  damping  ratio,  delta_h 
frequency  dependent  damping  ratio,  delta_i 
crack  length  dependent  damping  ratio,  delta_crack 
total  damping  ratio,  delta 
fatigue  loading,  p 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  1987 


double  *data; 
int  n; 
t 


/*  print  to  screen  */ 
printfl  "fatigue  cycle  nunber 
printf(  "fatigue  crack  length 
printfl  "stress  intensity  factor 
printfl  "torsional  stiffness  at  crack 
printfl  "fraction  of  life  expended 
printfl  “plastic  zone  size 
printfl  "plastic  zone  strain 
printfl 
printfl 
printfl 
printfl 
printfl 
printfl 


X5d  lcycle)\n",  n  ); 

X1 1.4 1 e  t in)\n",  datatO]  ); 

XII. 4le  tpsi - in‘1/2)\n",  datall]  ) 
XII. 4le  l lb- in)\n»,  datat2]  ); 
X11.4le\n",  data  13]  ); 

XII. 4le  lin)\n",  data[4]  ); 

XII. 4le  lin/in)\n“,  data[5]  ); 

XII. 4te  lrad/sec)\n",  data[6]  ); 
X11.4le\n",  data [7]  }; 

X11.4le\n",  data [8]  ); 

X11.4le\n“,  data [9]  >; 

X11.4le\n",  data  CIO]  ); 

XII. 4le  l lb)\n\n",  dataCII]  >; 


"resonant  frequency 
"amplitude-dependent  damping  ratio 
"frequency  dependent  damping  ratio 
"crack  length  dependent  damping  ratio 
"total  damping  ratio 
"fatigue  load 


return; 

) 


doctole  dddphit  x,  co  ) 

/*  a******************************************************************** 

NAME 

dddphi . third  derivative  of  beam  characteristic 

eauation  with  respect  to  x 

SYNOPSIS 


y  =  dddphi <  x,  co  ); 


double  *co, 

x. 

y; 


DESCRIPTION 


calculate  third  derivative  of  beam 
characteristic  equation  with 
respect  to  x 

coefficients  of  characteristic  equation 
beam  coordinate 

third  derivative  of  beam  at  location  x 


Calculate  third  derivative  of  beam  at  location  x.  Third  derivative 
has  the  following  equation  form 

dddphilx)  =  scale_factor*cotO] '3*1  sinl  co[0]*x  )  - 

co[1]*cosf  co[0J*x  )  ♦ 
co(2]*sinhl  co[0]*x  )  ♦ 
co[3]*coshf  co[0]*x  )  ) 


cotO]  =  kappa; 

coll]  =  mode. element [2]; 
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co [2]  =  mode. element [3] ; 

coC3]  =  mode. element  C4]  ; 
coU]  =  scale_factor; 

FUNCTIONS  USED 

dcos . lattice  SSP  cosine  function 

dcosh . Lattice  SSP  hyperbol  ic  cosine  function 

dsin . lattice  SSP  sine  function 

dsinh . Lattice  SSP  hyperbolic  sine  function 


REFERENCES 

Harris,  C.M.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  U.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 

*********************************************************************  * i 

double  *co,  x; 

{ 

double  y; 

double  dcosO,  dcoshO,  dsinO,  dsinhO; 

/*  calculate  third  derivative  of  beam  at  x  */ 

y  =  co[4]*cube(  co[0]  )*(  dsin(  co[0]*x  )  -  cot1]*dcos(  co[0l*x  )  + 
co[2]*dsinh(  co[0]*x  )  ♦  coC33*dcosh(  coCOI'x  }  ); 


return  (  y  ); 

> 


double  ddphi (  X,  co  ) 

j*  ********************************************************************* 

NAME 


ddr  ri . second  derivative  of  beam 

characteristic  equation  with  respect 
to  x. 


SYNOPSIS 

y  =  ddphi(  x,  co  );  calculate  second  derivative  of 

beam  characteristic  equation  with 
respect  to  x 

double  *co,  coefficients  of  characteristic  equation 

x,  beam  coordinate 

y;  curvature  of  beam  at  location  x 

DESCRIPTION 

Calculate  curvature  of  beam  at  location  x.  Curvature  has  the 
following  equation  form 

ddphi(x)  =  scale_factor*co[0] "2*(  -cos(  co[0]*x  )  - 

co[1]*sin(  co[0]*x  )  ♦ 
co[2]*cosh(  co[0]*x  )  ♦ 
co[3]*sinh(  cot0)*x  )  ) 

co[0]  =  kappa; 

cotlj  =  mode. element [2]; 
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co £2]  =  mode. element  [3]  ; 

coC3]  =  mode. element  [4]  ; 
co[4]  =  scale_factor; 

FUNCTIONS  USED 

dcos . Lattice  SSP  cosine  function 

dcosh . Lattice  SSP  hyperbolic  cosine  function 

dsin . Lattice  SSP  sine  function 

dsinh . Lattice  SSP  hyperbolic  sine  function 


REFERENCES 

Harris,  C.M.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  U.T.  "Theory  of  Vibration  with  Applications,"  Prentict- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 

A********************************************************************  *  j 

double  *co,  x; 

{ 

double  y; 

double  dcosO,  dcosh(),  dsin(),  dsinhO; 

/*  calculate  curvature  of  beam  at  x  */ 

y  =  co[4J*sqr(  co[0]  )*(  -dcos(  co[0]*x  )  -  co£1)*dsin(  co[0]*x  )  + 
co[2]*dcosh(  col0]*x  )  +  co[3]*dsinh(  co[0]*x  )  ); 


return  <  y  ); 

> 


dot±>le  ddphi2(  x,  co  ) 

/*  ********************************************************************* 
NAME 

ddphi2 . square  of  second  derivative  of  beam 

characteristic  equation  with  respect 
to  x. 


SYNOPSIS 

y  =  ddphi2(  x,  co  ); 


double  *co, 

x. 

y; 


DESCRIPTION 


calculate  square  of  second  derivative  of 
beam  characteristic  equation  with 
respect  to  x 

coefficients  of  characteristic  equation 
beam  coordinate 

curvature  of  beam  at  location  x 


Calculate  square  of  curvature  of  beam  at  location  x.  Curvature  has 
the  following  equation  form 

ddphi(x)  =  scale_factor*co[OJ ‘2*(  -cos(  col0]*x  )  - 

coC1J*sin(  col 0]*x  )  * 
co[2]*cosh(  co[0]*x  )  + 
co[3]*sinh(  co[0]*x  )  ) 


co[0]  =  kappa; 
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co til  s  mode. element [2] ; 
co[2]  =  mode. element [3] ; 

co 13]  =  mode. element [4] ; 

co[4]  =  scale_factor; 

FUNCTIONS  USED 

REFERENCES 

Harris,  C.H.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  U.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 

*************************** ******************************************  */ 


double  *co,  x; 

{ 

double  ddphiO; 

/*  calculate  square  of  curvature  of  beam  at  x  */ 
return  (  sqr(  ddphi(  x,  co  )  )  ); 

> 


doc±>le  dphi(  x,  co  ) 

NAME 

dphi . first  derivative  of  beam  characteristic 

equation  with  respect  to  x. 


SYNOPSIS 

y  =  dphi(  x,  co  );  calculate  first  derivative  of  beam 

characteristic  equation  with 
respect  to  x 

double  *co,  coefficients  of  characteristic  equation 

x,  beam  coordinate 

y;  slope  of  beam  at  location  x 

DESCRIPTION 

Calculate  slope  of  beam  at  location  x.  Curvature  has  the 

following  equation  form 

dphi(x)  =  scale_factor*colO] *(  -sin(  co[0]*x  )  + 

coC1]*cos(  co[0]*x  )  ♦ 
col2]*sinh(  co[0]*x  )  ♦ 
cot3]*cosh(  co[0j*x  )  ) 


co 10]  =  kappa; 

coCU  =  mode. element [2] ; 
co[2]  =  mode. element  13] ; 

co[3]  =  mode. element  14] ; 

co (4]  =  scale_factor; 

FUNCTIONS  USED 

dcos . Lattice  SSP  cosine  function 
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dcosh . Lattice  SSP  hyperbolic  cosine  function 

dsin . Lattice  SSP  sine  function 

dsinh . Lattice  SSP  hyperbolic  sine  function 


REFERENCES 

Harris,  C.M.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  U.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 

*********************************************************************  * j 

double  *co,  x; 

(. 

double  y; 

double  dcosO,  dcoshO,  dsinO,  dsinh{); 

/*  calculate  slope  of  beam  at  x  V 

y  =  co[4]*co[0J*(  -dsin(  colOI'x  )  +  coC1]*dcos<  co[0]*x  )  + 

co[2]*dsinh(  co[0]*x  )  +  co[3]*dcosh(  co[0]*x  )  ); 

return  (  y  ); 

> 


/•ARGSUSED  parameter  x  is  not  used  in  this  function  */ 
double  ezb_liarit(  x,  co  ) 

/*  ********************************************************************* 
NAME 

ezb_limit . z  coordinate  lower  limit  of 

integration  over  the  elastic  range 

SYNOPSIS 

y  =  ezb_limit(  x,  co  );  z  coordinate  lower  limit  of 

integration  over  the  elastic  beam 
material 

double  *co,  user  modifiable  coefficients  for 

defining  lower  limit 

x,  value  of  x  coordinate  for  application 

in  setting  lower  limit 

y;  lower  limit  of  integration 

DESCRIPTION 

This  function  defines  the  lower  limit  of  integration  for  the  z 
coordinate  of  a  beam.  The  limit  is  specifically  for  the  elastic 
range  of  material.  The  function  is  called  by  multi_al,  a  routine 
that  performs  multiple  numerical  integration. 

FUNCTIONS  USED 

REFERENCES 

Burden,  R.L.  and  Faires,  J.D.  "Nunerical  Analysis,"  third 
edition,  Prindle,  Weber  and  Schmidt,  Boston,  1984,  pp.  153- 
195. 


AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 

*********************************************************************  * / 

double  *co,  x; 

<. 

/*  the  first  element  in  the  array  of  coefficients  contains  a  constant 
that  defines  the  lower  limit  for  the  function  being  integrated  */ 

return  {  co[0]  ); 

> 


/‘ARGSUSED  parameter  x  is  not  used  in  this  function  */ 
double  eze_limit(  x,  co  ) 

/•  ***•****••..*****••**••*•••*•.**••*•*•.*.***••••••••••••**••*•••****• 

NAME 

eze_limit . z  coordinate  upper  limit  of 

integration  over  the  elastic  range 

SYNOPSIS 

y  =  ezb_limit(  x,  co  );  z  coordinate  upper  limit  of 

integration  over  the  elastic  beam 
material 

double  ‘co,  user  modifiable  coefficients  for 

defining  upper  limit 

x,  value  of  x  coordinate  for  application 

in  setting  upper  limit 

y;  upper  limit  of  integration 

DESCRIPTION 

This  function  defines  the  upper  limit  of  integration  for  the  z 
coordinate  of  a  beam.  The  limit  is  specifically  for  the  elastic 
range  of  material.  The  function  is  called  by  multi_al,  a  routine 
that  performs  multiple  numerical  integration. 

FUNCTIONS  USED 

REFERENCES 

Burden,  R.L.  and  Faires,  J.D.  "Numerical  Analysis,"  third 
edition,  Prindle,  Ueber  and  Schmidt,  Boston,  1984,  pp.  153- 
195. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 

*********************************************************************  • j 

double  *co,  x; 

< 

/*  the  second  element  in  the  array  of  coefficients  contains  a 
constant  that  defines  the  upper  limit  for  the  function  being 
integrated  */ 

return  (  coll]  ); 

) 
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int  get_data(  fp,  input_int,  input_real  ) 

; *  ********************************************************************* 
NAME 

get_data.... . input  data  from  file 

SYNOPSIS 

include  <stdio.h> 

err  =  get_data(  fp,  input_int,  input_real  );  input  material  and 

program  control  constants 

input  file 

pointer  to  array  of  real  constants 
error  code:  0  -  no  error 
1  -  error 

pointer  to  array  of  integer  constants 

DESCRIPTION 


FILE  *fp; 

double  *input_real; 
int  err, 

*input_int; 


get_data  accesses  an  open  file  and  inputs  the  data  associated  with 
the  two  arrays.  Data  is  placed  in  the  input  file  in  the  following 
format: 


input_real : 

a_0 
a_disl 
a_f i na l 
alpha 
b_crack 

c 

c  1 
c~2 
c_heat 
chi 

delta_ext 

e 

ei 

n_max 

n_ult 

g_check 

growth_rate 


h 

k 

k_kni fe_edge 
l 

l_exciter 

l_load 

lambda_0 

m_beam 

m_exciter 

p” 


=  initial  crack  length  (in). 

=  dislocation  density  coefficient  (disl/in‘2). 

=  final  crack  length  (in). 

=  linear  expansion  coefficient  (in/in/xF). 

=  material  constant  relating  crack  length  to 
torsional  stiffness  (lb- in). 

=  Manson-Coff in  ductility  exponent. 

=  internal  damping  coefficient. 

=  internal  damping  coefficient  (in/in). 

=  specific  heat  per  unit  volume  (psi/xF) 

=  strain  independent  probability  for  the 
remobilization  of  immobile  dislocations. 

=  damping  ratio  caused  by  external  source 
=  Young's  modulus  (psi). 

=  Young's  modulus  x  moment  of  inertia  (lb-in‘2). 
=  maximum  strain  in  beam  during  damping  data 
collection  (in/in). 

=  static  plastic  fracture  value  at  1/2  cycle 
(in/in). 

=  minimum  crack  growth  before  updating 
eigenvector  results  (in). 

=  crack  growth  rate  constant  based  on 
dislocation  theory  of  crack  growth 
(in'7/lb‘4-cycle). 

=  beam,  height  (in). 

=  thermal  conductivity  (Ib/sec-xF) 

=  torsional  stiffness  at  knife-edge  support 
(lb-in). 

=  beam  length 

=  location  of  exciter  relative  to  knife-edge 
support  (in). 

=  location  of  load  relative  to  knife-edge 
support  (in). 

=  dislocatin  density  for  n_p  =  0  (disl/in‘2). 

=  mass  of  beam  per  unit  length  (slug/inch). 

=  mass  of  exciter  and  counterweight  (slug). 

=  fatigue  load  applied  to  beam  (lb). 
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e_yield 

slope_factor 


=  yield  stress  (psi). 

=  crack  dependent  constant  relating  beam  slope 
and  damping. 

temperature  =  ambient  temperature  of  specimen  (xF). 
toler  =  tolerance  for  iteration  algorithms. 

mode. element [0]  thru  £41  are  initial  guesses  for  the  newton 
iterative  solution  to  the  nonlinear  beam  equation,  and 
represent  the  following  data 


elementlO]  -  eigenvalue,  kappa*l 
element £1]  -  resonant  frequency,  omega 
element £2]  -  B  coefficient 
element £3]  -  C  coefficient 
element £4]  -  D  coefficient 

input_int: 

data_rate  =  nunber  of  fatigue  cycles  between  data  output 
is_type_sif  =  0  for  constant  load 

1  for  constant  stress  intensity  factor 
j  =  initial  cycle  nurber 

m  =  number  of  subintervals  for  integration 

maxum  =  number  of  iteration  cycles 

n  =  number  of  subintervals  for  second  integration 

FUNCTIONS  USED 

fscanf . read  from  file 


REFERENCES 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  19  AUG  87 

at********************************************************************  »j 
FILE  *fp; 

double  *input_real; 
int  *input_int; 
t 

int  err,  i; 

err  =  0;  /*  initialize  return  code  */ 

/*  input  all  real  values,  one  per  line  */ 
for  (  i  =  0;  i  <  MAX  REAL;  i++  ) 

{ 

if  (  fscanf(  fp,  "%lf",  (  input_real  +  i  )  )  ==  0  ) 

{ 

/*  error  in  input  */ 

printf(  "File  error  reading  reals  in  get_data,  i  =  Xd\n",  i  ); 
err  =1;  /*  set  error  code  V 

break;  /*  end  for  loop  */ 

> 

> 

/*  input  all  integer  values,  one  per  line  */ 
if  (  err  ==  0  ) 

< 

for  (  i  =  0;  i  <  MAX  INT;  i++  ) 

{ 

if  (  fscanf(  fp,  "%d",  (  input  int  ♦  i  )  )  ==  0  ) 

< 

/*  error  in  input  */ 

printf(  "File  error  reading  integers  in  get_data,  i  =  Xd\n",  i  ) 
err  =1;  /*  set  error  code  */ 

break;  /*  end  for  loop  */ 
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> 

> 

) 

/*  return  error  code  V 
return  (  err  ); 


void  input_print(  input  int,  input  real  ) 

/*  a*****....*****.****;******..**;...**..******.******.****.***...*.*.. 

NAME 

inputjarint . print  input  data  to  screen 

SYNOPSIS 

include  <stdio.h> 

input_print(  input_int,  input_real  );  print  input  material  and 

program  control  constants  to  screen 

double  *input_real;  pointer  to  array  of  real  constants 

int  *input_int;  pointer  to  array  of  integer  constants 

DESCRIPTION 


This  function  takes  the  two  input  arrays  and  prints  them  on  the 
the  screen  along  with  their  descriptions. 


input_real : 

a_0 
a_disl 
a_f inal 
alpha 
b_crack 

c 

C_1 

c_2 

c_heat 

chi 

delta_ext 

e 

ei 

njnax 

n_ult 

g_check 

growth_rate 


h 

k 

k_kni  fe__edge 
l 

l_exciter 

l_load 

lambda_0 
m_beam 
m  exciter 


=  initial  crack  length  (in). 

=  dislocation  density  coefficient  (disl/in‘2). 

=  final  crack  length  (in). 

=  linear  expansion  coefficient  (in/in/xF). 

=  material  constant  relating  crack  length  to 
torsional  stiffness  (lb- in). 

=  Manson-Coff in  ductility  exponent. 

=  internal  damping  coefficient. 

=  internal  damping  coefficient  (in/in). 

=  specific  heat  per  unit  volume  (psi/xF) 

=  strain  independent  probability  for  the 
remobilization  of  immobile  dislocations. 

=  damping  ratio  caused  by  external  source 
=  Young's  modulus  (psi). 

=  Young's  modulus  x  moment  of  inertia  (lb-in'2). 
=  maximum  strain  in  beam  during  damping  data 
collection  (in/in). 

=  static  plastic  fracture  value  at  1/2  cycle 
(in/in). 

=  minimum  crack  growth  before  updating 
eigenvector  results  (in). 

=  crack  growth  rate  constant  based  on 
dislocation  theory  of  crack  growth 
(in‘7/lb'4-cycle) . 

=  beam  height  (in). 

=  thermal  conductivity  (Ib/sec-xF) 

=  torsional  stiffness  at  knife-edge  support 
(lb-in). 

=  beam  length 

=  location  of  exciter  relative  to  knife-edge 
support  (in). 

=  location  of  load  relative  to  knife-edge 
support  (in). 

=  dislocatin  density  for  n_p  =  0  (disl/in‘2). 

=  mass  of  beam  per  unit  length  (slug/inch). 

=  mass  of  exciter  and  counterweight  (slug). 
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p 

e_yield 

sTope_factor 

temperature 

toler 


=  fatigue  toad  applied  to  beam  (lb). 

=  yield  stress  (psi). 

=  crack  dependent  constant  relating  beam  slope 
and  damping. 

=  ambient  temperature  of  specimen  (xF). 

=  tolerance  for  iteration  algorithms. 


mode. element  CO)  thru  [4]  are  initial  guesses  for  the  newton 
iterative  solution  to  the  nonlinear  beam  equation,  and 
represent  the  following  data 


element [01  -  eigenvalue,  kappa*l 
element [1]  -  resonant  frequency,  omega 
elemental  -  8  coefficient 
element [3]  -  C  coefficient 
elemental  -  0  coefficient 


input_int: 

data_rate 

is_type_sif 

j 

m 

maxim 

n 

FUNCTIONS  USED 
REFERENCES 


=  number  of  fatigue  cycles  between  data  output 
=  0  for  constant  load 

1  for  constant  stress  intensity  factor 
=  initial  cycle  nuiber 
=  number  of  subintervals  for  integration 
=  nuiber  of  iteration  cycles 

=  number  of  subintervals  for  second  integration 


AUTHOR:  Rick  Jones 

VERSION:  1.0 
OATE:  19  AUG  87 


V 


double  *input_real; 
int  "input  int; 

( 

/*  print  real  input  data  to  screen  */ 

printf(  "initial  crack  length 

printf(  “dislocation  density  coefficient 

printf(  "final  crack  length 

printf(  "linear  expansion  coefficient 

printf(  "crack  torsional  stiffness  slope 

printf(  "Nanson-Coff in  ductility  exponent 

printf(  "internal  damping  coefficent  1 

printf(  "internal  damping  coefficent  2 

printff  "specific  heat 

printf(  "remobilization  probability 

printf(  "external  damping  ratio 

printf(  "Young's  modulus 

printf(  "modulus  x  Moment  of  Inertia 

printfC  "damping  cycle  maximum  strain 

printf(  "static  plastic  fracture  value 

printf(  "eigenvalue  update  parameter 

printf(  "crack  growth  rate 

printf(  "beam  height 

printf(  "thermal  conductivity 

printf(  "knife-edge  torsional  stiffness 

printf(  "beam  length 

printf(  "excitation  location 

printf(  "initial  dislocation  density 

printfC  “fatigue  load  location 

printf(  "beam  mass  per  unit  length 

printf(  "exciter  mass 

printff  "fatigue  load 

printf(  "yield  stress 


XII. 4le  (in)\n",  input_real [0]  ); 

XII. 4le  (disl/in'2)\nH,  input_real Cl]  ); 

XII. 4le  (in)\n",  input_real [2]  ); 

XII. 4le  (in/in/xF)\n",  input_real C3]  ); 

XII. 4le  (lb-in)\nH,  input_real [4]  ); 

X11.4le  \n",  input_real E5T  ); 

XII. 4le  \n",  input_real C6]  ); 

XII. 4le  (in/in)\n",  input_real [7]  ); 

XII. 4le  (psi/xF)\n",  input_real C8]  ); 
X11.4le  \n",  input_real C9]  ); 

XII. 4le  \n",  input_real CIO]  ); 

XII. 4le  (psi)\n",  input_real [11]  ); 

XII. 4le  ( lb- in‘2)\n",  input_real [12]  ); 

XII. 4te  (in/in)\n",  input_real [13J  ); 

XII. 4le  (in/in)\n",  input_real [14]  ); 

XII. 4le  Cin)\n",  input_real [15]  ); 

XII. 4le  (in'7/lb‘4-cycTe)\n",  input_real [16] 
XII. 4le  (in)\n",  input_real [17]  ); 

XII. 4le  <lb/sec-xF)\n",  input_real [18]  >; 
XII. 4le  <lb-in)\n",  input_real [19]  ); 

XII. 4le  (in)\n",  input_real [20]  ); 

XII .4le  ( in)\n",  input_real [21]  ); 

XII. 4le  (disl/in'2)\n",  input_real [22]  ); 
XII. 4le  (in)\n",  input_real [23]  ); 

XII. 4le  (slug/in)\n",  input_real [24]  ); 

XII. 4le  (slug)\n",  input_real [25]  ); 

XII. 4le  (lb)\n",  input_real [26]  ); 

XII. 4le  (psi)\n",  input_real [27]  ); 
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printf(  "damping  slope  factor 
printf(  "temperature 
printff  "approximation  tolerance 
printff  "initial  guess  of  beam  mode 
printff  “  eigenvalue 

printf(  "  resonant  frequency 

printfl  "  B  coefficient 

printf(  "  C  coefficient 

printfl  "  D  coefficient 


X11.4te  \n",  input_real [28]  ); 

XII. 4le  (xF)\n",  input_real 129]  ); 
X11.4le  \n",  input_real[30]  ); 
parameters\n"  ); 

XII. 4le  \n",  input_real  t31]  ); 

X11.4le  (rad/sec)\n",  input_real 132]  ); 
X1 1.4 1 e  \n",  input_real [33]  ); 

XII. 4le  \n",  input_real [34]  ); 

XII. 4le  \n",  input_real [35]  ); 


/*  print  integer  input  data  to  screen  */ 
printf(  "sampling  interval  X5d 


(cycle)\n",  input_int[0]  ); 


/*  report  load  type  */ 
if  (  *(  input_int  +  1  )  ==  0  ) 

printf(  “load  type  -  constant  load\n"  ); 
else 

printff  "load  type  -  constant  stress  intensity  factor\n"  ); 
printf(  "initial  cycle  nurber  X5d 

printf(  "subintervals  for  integration  X5d 

printf(  "maximun  iteration  cycles  X5d 


printft  "subintervals  on  2nd  integration  XSd 
return; 


(cycle)\n",  input_int[2]  ); 
(intervals)\n",  input_int[3] 
<cycle)\n",  input_intT4]  ); 


>; 


(intervals)\n",  input_int[5]  ); 


docfcle  f(  r  ) 

!*  tit***************************** ************************************** 

NAME 

f . empirical  formula  for  pure  bending 

stress  intensity  factor 

SYNOPSIS 

double  f(); 

y  =  f(  r  >;  pure  bending  stress  intensity  factor 

empirical  formula 

double  r,  ratio  of  crack  length  to  section  depth 

y;  empirical  formula  for  pure  bending 

DESCRIPTION 

Calculate  the  empirical  formula  for  pure  bending 

f(a/h)  =  1.122  -  1.4*(a/h)  ♦  7.33*(a/h)'2  •  13.08*<a/h)‘3  + 
14*(a/h)‘4 

(  0.2X  for  a/h  <  0.6), 
or 


f (a/h)  =  sqrt(tan(p)/p)* 

(.923  ♦  . 199*( 1  -  sin(p)))'4/cos(p), 

where  p  =  c*a/(2*h)  and  r  =  a/h 

(  0.5X  for  any  a/h). 

FUNCTIONS  USED 


dcos . lattice  SSP  cosine  function 

dsin . lattice  SSP  sine  function 

dsqrt . Lattice  SSP  square  root  function 

dtan . lattice  SSP  tangent  function 
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REFERENCES 


Tada,  M.,  Paris,  P.C.  and  Irwin,  G.R.,  "The  Stress  Analysis  of 
Cracks  Handbook,"  Del  Research  Corporation,  Hellertown, 
Pennsylvania,  1973,  pp  2. 13-2. 14. 


AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  19  AUG  87 


*/ 


double  r; 

( 

double  p; 

double  dcosO,  dsinO,  dsqrtO,  dtanO; 

/*  select  the  formula  */ 
if  C  r  <=  0.6  ) 

/*  return  polynomial  calculation  */ 

return  (  1.122  -  1.4*r  ♦  7.33*sqr<  r  )  -  13.08*cube(  r  )  ♦  14.*quad(  r  )  ); 
else  { 

/*  return  trigonometric  calculation  */ 
p  =  .5*PI*r; 

return  (  dsqrt(  dtan(  p  )/p  )*quad(  .923  ♦  . 199*(  1.  -  dsin(  p  )  )  )/dcos(  p  )  ) 

> 


double  ki(  a,  h,  e  ) 

/*  *************•*******•***•*•****.••**••.*«**.«*******••••.••**»•*•. 
NAME 

ki . stress  intensity  factor  for  pure 

bending 

SYNOPSIS 

y  =  ki <  a,  h,  e  );  stress  intensity  factor  for  the  case  of 

pure  bending 

crack  length 
section  depth 
nominal  stress 
stress  intensity  factor 

DESCRIPTION 


double  a, 
h. 
e. 
y. 


Begin  by  calculating  the  stress  intensity  factor  using  an 
empirical  formula  for  pure  bending 

f(e/h)  =  1.122  -  1 ,4*(a/h)  ♦  7.33*<a/h)‘2  -  13.08*Ca/h)‘3  ♦ 
14*(a/h)'4 

(  0.2X  for  a/h  <  0.6), 
or 

f(a/h)  =  sqrt(tan(psi )/psi )* 

(.923  +  .199*0  -  sin(psi)))"4/cos(psi), 

where  psi  =  c*a/(2*h) 

(  0.5X  for  any  a/h). 

y  =  e*sqrt(c*a)*f(a/h)  stress  intensity  factor 
FUNCTIONS  USED 
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dsqrt.. 

REFERENCES 


Lattice  SSP  square  root  function 


Tada,  H.,  Paris,  P.C.  and  Irwin,  G.R.,  "The  Stress  Analysis  of 
Cracks  Handbook,"  Del  Research  Corporation,  Hetlertown, 
Pennsylvania,  1973,  pp  2.13-2.14. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  19  AUG  87 

.**.**.**...••*.*****...**..*****.*****••***.***..**************•*•.. 

double  a,  h,  e; 

< 

double  f<),  dsqrt (); 

/*  return  the  stress  intensity  factor  */ 
return  (  e*dsqrt(  PI*a  )*f(  a/h  )  ); 


MATRIX  Bodal(  x,  co,  residuals  ) 

j*  ********************************************************************* 

NAME 

modal . calculates  residuals 

SYNOPSIS 

^include  <numeric.h> 

residuals  =  modal(  x,  co,  residuals  );  calculate  residuals  of 

nonlinear  equations 

double  *co;  pointer  to  array  of  nonlinear  equation 

coefficients 

MATRIX  residuals  vector  containing  equation  errors  based 

on  equation  coefficients  and  vector 
of  guess  values 

x;  vector  of  guess  values  to  solve 

nonlinear  set  of  equations 

DESCRIPTION 

Contains  user  supplied  equations  to  determine  residuals  in  solving 
nonlinear  system  of  equations  in  x.  This  function  is  designed  for 
use  with  newtonO,  the  Newton-Raphson  Iterative  process  for 
solving  nonlinear  sets  of  equations. 


cotO]  =  k_knife_edge; 
coll]  =  k_crack; 
co[2]  =  l; 
co  13]  =  m_beam; 
co [4]  =  m_exciter; 
cot5]  =  el; 
col6]  =  l_exciter; 

cot7]  =  (double)  m; 


knife-edge  stiffness 
crack  stiffness 
beam  length 

beam  mass  per  unit  length 

mass  of  exciter 

modulus  and  moment  of  inertia 

position  of  exciter  relative  to  knife- 

edge 

number  of  subdivisions  for  integration 


FUNCTIONS  USED 


coeff_red . reduce  characteristic  oefficients 

dcos . Lattice  SSP  cosine  function 


dcosh . Lattice  SSP  hyperbolic  cosine  function 

ddphi2 . square  of  the  second  derivative  of  the 

characteristic  equation 

dsin . Lattice  SSP  sine  function 

dsinh . Lattice  SSP  hyperbolic  sine  function 

mat_assign . assign  matrix  parameters  to  variable 

mat_delete . delete  matrix  storage  space 

phi2 . square  of  the  characteristic  equation 

simpson . numerical  integration  using  Simpson's  rule 


*/ 


double  *co; 

MATRIX  residuals,  x; 

< 

D_SUPPLD  h_fun,  ddh_fun; 

double  ch_ei,  c_ei,  eigen,  gen_mass,  gen_stiff; 

double  sh_ei,  s_ei,  stiff_c,  stiff_ke; 

double  phi_data75]; 

double  ddphi2(),  phi2(); 

double  dcosO,  dcoshO,  dsin(),  dsinhO; 

int  i; 

MATRIX  a,  coeff; 

MATRIX  coeff_red(),  mat_assign(),  mat_detete(); 
/*  user  supplied  equations  for  beam  analysis  */ 


eigen 

ch_ei 

c_ei 

sh_ei 

s_ei 

stiff_c 

stiff  ke 


=  x.element!0]*co!2];  /*  eigenvalue  *  length  */ 

=  dcosh (  eigen  ); 

=  dcos!  eigen  ); 

=  dsinh(  eigen  ); 

=  dsin(  eigen  ); 

=  co[1]/(  x. element I0]*co!5]  );  /*  characteristic  stiffness 

at  crack  */ 

=  colO]/(  x. element [0]*col5]  );  /*  characteristic  stiffness 

at  boundary  */ 


hfun.funct  =  phi2; 

hfun. coeff  =  phi_data; 


/*  beam  characteristic  equation 
/*  eigenvector  data 


V 

V 


ddh_fun.funct  =  ddphi2; 
dcfti_fun. coeff  =  phi_data; 


/*  second  derivative  of  beam  characteristic 
equation  */ 

/*  eigenvector  data  */ 


/*  characteristic  equation  coefficients  */ 
phi_data[0]  =  x. element [0] ; 
phi_data[1]  =  x. element  [2J  ; 
phi_data[2]  =  x. element  13] ; 
phi_data[3]  =  x.  element  14]  ; 

phi_datal4]  =1.;  /*  use  scale  factor  of  1  */ 


/*  create  a  matrix  */ 
a  =  mat_assign(  4,  4,  NULL  ); 


/*  fill  a  matr 
a. element [  0] 
a. element!  1] 
a. element!  2] 
a. element!  3] 
a. element!  4] 
a. element!  5] 
a. element!  6] 
a. element!  71 
a. element!  8] 
a. element !  9] 
a. element  110] 
a. element !11] 
a. element  112] 
a. element !13] 


ix  */ 

=  i.; 

=  1.; 

=  -(  c_ei  +  stiff_c*s_ei 
=  s_ei ; 

=  0.; 

=  stiff_ke; 

=  -s_ei  ♦  stiff_c*c_ei; 

=  -c  ei; 

■  1.7 

=  -1.; 

=  ch_ei  +  st i f t_c*sh_ei ; 
=  sh_ei; 

=  0.; 

=  stiff_ke; 


); 
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a. element [141  =  sh_ei  +  stiff_c*ch_ei; 
a. element [IS]  =  ch_ei; 

/*  determinate  of  a  should  be  0,  first  residual  is  det(a)  */ 
res i duals. el ement[0]  =  mat_determ<  a  ); 

/*  create  a  vector  to  hold  coefficients  */ 
coeff  *  mat_assign(  4,  1,  NULL  ); 

/*  fill  coefficients  in  terms  of  coeff [0]  */ 
coeff  =  coeff_red(  a,  coeff  ); 

/*  calculate  residuals  for  mode  shape  coefficents  */ 
for  (  i  *  1;  i  <*  3;  ++i  ) 

residuals. element [i*1]  *  coeff .element [i]  -  x. element [i+1] ; 

/*  beam  mass  */ 

genjnass  =  co[3]*simpson(  0.,  co[2],  (int)  co[71,  h_fun  )  + 
co[4]*phi2(  co[6],  phi_data  ); 


/*  beam  stiffness  */ 

gen_stiff  *  co(5]*simpson(  0.,  co[21,  (int)  co[7],  ddh_fun  ); 
/*  residual  for  frequency  */ 

residuals. element [1]  =  sqr(  x.element[1]  )*gen_mass  -  gen_stiff; 

/*  release  storage  */ 
a  =  mat_delete(  a  ); 

coeff  =  mat_delete(  coeff  ); 

return(  residuals  ); 

) 


docfcle  n2_theta_h(  x,  z,  co  ) 
NAME 


n2_theta_h . beam  strain  squared  multiplied  by 

internal  damping  ratio 

SYNOPSIS 

y  =  n2_theta_h(  x,  z,  co  );  beam  strain  squared  times  internal 

damping  ratio 

expression  coefficients 
length  coordinate  of  beam 
result 

height  coordinate  of  beam 

DESCRIPTION 


double  *cof 

y. 

z; 


Calculate  beam  strain  squared  multiplied  by  internal  damping  ratio 
at  any  x,  z  coordinate  of  a  beam.  The  result  is  computed  from  the 
following  equation 

y  =  n(x,z)‘2*c_1*(c_2/n(x,z))'2*ei1(c_2/n<x,z>), 

where  c_1  and  c_?  are  material  constants.  n()  is  the  beam  strain 
and  ei1()  is  the  exponential  integral  function. 


co[0] 

=  O; 

internal  damping  coefficient 

com 

=  c_2; 

internal  damping  coefficient 

co  [2] 

=  kappa; 

eigenvalue 

co[3] 

=  mode. element [21 ; 

phi(x)  coefficient  B 

co  [41 

=  mode. element [31; 

phi(x)  coefficient  C 

co[51 

=  mode. element [4]; 

phi(x)  coefficient  D 
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co  [6]  =  phi_data[4]; 


amplitude  scale  factor 


FUNCTIONS  USED 


eil . Lattice  SSP  exponential  integral 

function 

strain . beam  strain 


REFERENCES 

Harris,  C.M.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  U.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

Povolo,  F.  "On  the  Granato-Lucke  Expression  for  the  Amplitude- 
Dependent  Damping,"  Scripta  Metallurgies,  v  9,  pp  865-872,  1975. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  23  AUG  87 

*********************************************************************  * i 

double  *co,  x,  z; 

{ 

double  ratio; 

double  eilO,  strainO; 

int  ire; 

/*  initialize  error  code  */ 
ire  =  0; 

/*  test  input  to  exponential  integral  function  */ 
if  (  z  )  /*  z  coordinate  must  be  nonzero  */ 

ratio  =  abs(  co[1]/strain(  x,  z,  &co[2]  )  ); 
else 

return(  0.  ); 

if  <  ratio  >  669.  )  /*  test  for  underflow  */ 

return  (  0.  ); 
else 

/*  return  result  -  ignore  first  two  coefficients  when  providing 
pointer  to  strain  function  */ 
returnl  co[0]*sqr(  coll]  )*ei1(  ratio,  &irc  )  ); 


double  phi(  x,  co  ) 

I *  ********************************************************************* 

NAME 


phi . beam  characteristic  equation 

SYNOPSIS 

y  =  phi(  x,  co  );  calculate  displacement  from 

beam  characteristic  equation 

double  *co,  coefficients  of  characteristic  equation 

x,  beam  coordinate 

y;  displacement  of  beam  at  location  x 


DESCRIPTION 
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Calculate  displacement  of  beam  at  location  x.  Displacement  has 
the  following  equation  form 

phi(x)  =  scale_factor*(  cost  co[03*x  )  ♦  co[13*sin(  co[03*x  )  + 
co[23*cosh<  co[03*x  )  ♦  co[33*sinh(  co[03*x  )  ) 


co [03  =  kappa; 

co  Cl  3  =  mode. element [2] ; 

co[23  *  mode. element [33; 
co[33  =  mode. element  [43 : 

co[43  =  scale_factor; 


FUNCTIONS  USED 


dcos . Lattice  SSP  cosine  function 

dcosh . Lattice  SSP  hyperbolic  cosine  function 

dsin . Lattice  SSP  sine  function 

dsinh . Lattice  SSP  hyperbolic  sine  function 


REFERENCES 

Harris,  C.N.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  W.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 


V 


double  *co,  x; 

< 

double  y; 

double  dcos(),  dcoshO,  dsin(),  dsinh(); 

/*  calculate  displacement  of  beam  at  x  */ 
y  *  co[41*<  dcos(  co(03*x  )  ♦  co[1]*dsin(  co[03*x  )  ♦ 

co[2]*dcosh(  co[03*x  )  ♦  co[33*dsinh{  co[0]*x  )  ); 


return  (  y  ); 

> 


dodble  pfii2(  x,  co  ) 
/*  *****•***••*•••*• 
NAME 


phi2 . square  of  beam  characteristic  equation 

SYNOPSIS 

y  =  phi2(  x,  co  );  calculate  square  of  displacement  from 

beam  characteristic  et^tion 

double  *co,  coefficients  of  characteristic  equation 

x,  beam  coordinate 

y;  displacement  of  beam  at  location  x 


DESCRIPTION 

Calculate  square  of  displacement  of  beam  at  location  x.  Displacement 
has  the  following  equation  form 

phi(x)  =  scale_factor*(  cos(  co[0]*x  )  ♦  co[13*sin(  co[0]*x  )  ♦ 
co[23*cosh<  co[0]*x  )  ♦  co[33*sirV?(  col 03 *x  )  ) 
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co CO]  =  kappa; 
cod]  =  mode,  element  [2]; 
cot2]  =  mode. element [33 ; 
co[31  =  mode. element  [4]; 

co[ 4]  =  scale_factor; 

FUNCTIONS  USED 

phi . beam  characteristic  equation 

REFERENCES 

Harris,  C.M.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  W.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 


double  *co,  x; 

< 

double  phi(); 

/*  calculate  square  of  displacement  of  beam  at  x  */ 
return  (  sqr(  phi(  x,  co  )  )  ); 

> 


/•ARGSUSED  parameter  x  is  not  used  in  this  function  */ 
ckxijle  pzb  l  i*i  t(  x,  co  ) 

/*  ********************************************************************* 
NAME 

pzb_limit . z  coordinate  lower  limit  of 

integration  over  the  plastic  range 

SYNOPSIS 

y  =  pzb_limit(  x,  co  );  z  coordinate  lower  limit  of 

integration  over  the  plastic  beam 
material 

double  *co,  user  modifiable  coefficients  for 

defining  lower  limit 

x,  value  of  x  coordinate  for  application 

in  setting  lower  limit 

y;  lower  limit  of  integration 

DESCRIPTION 

This  function  defines  the  lower  limit  of  integration  for  the  z 
coordinate  of  a  beam.  The  limit  is  specifically  for  the  plastic 
range  of  material.  The  function  is  called  by  mutti_al,  a  routine 
that  performs  multiple  nunerical  integration. 

FUNCTIONS  USEO 

REFERENCES 
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Burden,  R.L.  and  Faires,  J.O.  "Numerical  Analysis,"  third 
edition,  Prindle,  Weber  and  Schmidt,  Boston,  1984,  pp.  153- 
195. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 

*********************************************************************  * ^ 

double  *co,  x; 

( 

/*  the  first  element  in  the  array  of  coefficients  contains  a  constant 
that  defines  the  lower  limit  for  the  function  being  integrated  */ 

return  (  co[0]  ); 

> 


/*ARGSUSED  parameter  x  is  not  used  in  this  function  */ 
docfcle  pzel ini t(  x,  co  ) 

NAME 

pze_limit . z  coordinate  upper  limit  of 

integration  over  the  plastic  range 

SYNOPSIS 

y  =  pzb_limit(  x,  co  );  z  coordinate  upper  limit  of 

integration  over  the  plastic  beam 
material 

double  *co,  user  modifiable  coefficients  for 

defining  upper  limit 

x.  value  of  x  coordinate  for  application 

in  setting  upper  limit 

y;  upper  limit  of  integration 

DESCRIPTION 

This  function  defines  the  upper  limit  of  integration  for  the  z 
coordinate  of  a  beam.  The  limit  is  specifically  for  the  plastic 
range  of  material.  The  function  is  called  by  multi_al,  a  routine 
that  performs  multiple  nunerical  integration. 

FUNCTIONS  USED 

REFERENCES 

Burden,  R.L.  and  Faires,  J.D.  "Numerical  Analysis,"  third 
edition,  Prindle,  Weber  and  Schmidt,  Boston,  1984,  pp.  1a3- 
195. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  18  AUG  87 

************* ***************** ***************************************  * y 

double  *co,  x; 

/*  the  second  element  in  the  array  of  coefficients  contains  a 
constant  that  defines  the  upper  limit  for  the  function  being 
integrated  */ 


D-42 


return  (  co£1]  ); 


> 


j *  ********************************************************************* 

NAME 

p2?? . partial  derivative  functions 

SYNOPSIS 

y  =  p2??(  x,  co  );  partial  deriviative  of  second  equation 

double  *co,  pointer  to  array  of  coefficients 

x,  location  on  beam 

y;  partial  derivative 

DESCRIPTION 

These  functions  are  partial  derivatives  of  beam  generalize  mass 
and  stiffness  with  respect  to  the  eigenvalue  and  the  beam 
characteristic  equation  coefficients: 

p2k . partial  of  generalized  stiffness  with 

respect  to  the  eigenvalue 

p2kb . partial  of  generalized  stiffness  with 

respect  to  the  B  coefficient 

p2kc . partial  of  generalized  stiffness  with 

respect  to  the  C  coefficient 

p2kd . partial  of  generalized  stiffness  with 

respect  to  the  D  coefficient 

p2m . partial  of  generalized  mass  with 

respect  to  the  eigenvalue 

p2mb . partial  of  generalized  mass  with 

respect  to  the  B  coefficient 

p2mc . partial  of  generalized  mass  with 

respect  to  the  C  coefficient 

p2md . partial  of  generalized  mass  with 

respect  to  the  D  coefficient 

co[0]  =  kappa; 

coll]  =  mode. element [2]; 

co[2]  =  mode. element C3] ; 

co 13]  =  mode. element [4] ; 

co[4]  =  scale_factor; 

FUNCTIONS  USED 

dcos . Lattice  SSP  cosine  function 

dcosh . Lattice  SSP  hyperbolic  cosine  function 

dsin . Lattice  SSP  sine  function 

dsinh . Lattice  SSP  hyperbolic  sine  function 

dphi . derivative  of  characteristic  equation 

with  respect  to  length  coordinate 

ddphi . second  derivative  of  characteristic 

equation  with  respect  to  length 
coordinate 

dddphi . third  derivative  of  characteristic 

equation  with  respect  to  length 
coordinate 

Phi .  . beam  characteristic  equation 

REFERENCES 

Burden,  R.L.  and  Faires,  J.D.  "Nunerical  Analysis,"  third 
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edition,  Prindle,  Weber  and  Schmidt,  8oston,  1984,  pp.  153- 
195. 


Harris,  C.M.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 


Thomson,  W.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 

AUTHOR:  Rick  Jones 

VERSION:  1.0 
DATE:  25  AUG  87 


a********************************************************************  */ 


double  p2k(  x,  co) 
double  *co,  x; 

{ 

double  ddphiO,  dddphiO; 

/*  partial  derivative  of  generalized  stiffness  with  respect  to 
eigenvalue  */ 

returnl  2.*x*ddphi(  x,  co  )*dddphi(  x,  co  )/co[0]  ); 


> 


double  p2kb(  x,  co) 
double  *co,  x; 

< 

double  ddphiO,  dsin(); 

/*  partial  derivative  of  generalized  stiffness  with  respect  to 
coefficient  B  */ 

return  (  -2-*sqr(  co[0]  )*ddphi(  x,  co  )*dsin(  col0]*x  )  ); 


> 


double  p2kc(  ,r,  co) 
double  *co,  x; 

{ 

double  ddphiO,  dcosh(); 

/*  partial  derivative  of  generalized  stiffness  with  respect  to 
coefficient  C  */ 

return  (  2.*sqr(  co[0]  )*ddphi(  x,  co  )*dcosh(  co[0]*x  )  ); 


> 


double  p2kd(  x,  co) 
double  *co,  x; 

double  ddphiO,  dsinhO; 

/*  partial  derivative  of  generalized  stiffness  with  respect  to 
coefficient  D  */ 

return  (  2.*sqr(  co(0]  )*ddphi(  x,  co  )*dsinh<  col0]*x  )  ); 


> 
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double  p2a(  x,  co) 
double  *co,  x; 

{ 

double  dphi(),  phi(); 

/*  partial  derivative  of  generalized  mass  with  respect  to  eigenvalue  */ 
return  (  2.*x*phi(  x,  co  )*dphi(  x,  co  )/co[0]  ); 


double  p2*b(  x,  co) 
double  *co,  x; 

{ 

double  dsinO,  phi(); 

/*  partial  derivative  of  generalized  mass  with  respect  to 
coefficient  S  */ 

return  (  2.*phi(  x,  co  )*dsin<  coI0]*x  )  ); 


> 


double  p2ac(  x,  co) 
double  *co,  x; 

< 

double  dcoshO,  pb it); 

/*  partial  derivative  of  generalized  mass  with  respect  to 
coefficient  C  */ 

return  (  2.*phi(  x,  co  )*dcosh(  co[0]*x  )  ); 


> 


double  p2ad(  x,  co) 
double  *co,  x; 

C 

double  dsinhO,  phiO; 

/*  partial  derivative  of  generalized  mass  with  respect  to 
coefficient  D  */ 

return  (  2.*phi(  x,  co  )*dsinh(  co[0]*x  )  ); 


> 


double  strainf  x,  z,  co  ) 
/*  ********************** 
NAME 


strain . beam  strain 

SYNOPSIS 


y  =  strainf  x,  z,  co  ); 
double  *co, 

x. 

Y. 

*; 


compute  beam  strain 

coefficients  used  in  strain  term 
x  coordinate  of  beam 
beam  strain 
z  coordinate  of  beam 


DESCRIPTION 
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Compute  strain  at  beam  coordinates  given  at  x  and  z. 
The  strain  is  computed  from  the  following  formula 


n(x,z)  =  z*phi"(x) 

co 10]  =  kappa; 

coll]  =  mode. element [2] ; 
co[2]  *  mode. element  13] ; 

cot3]  =  mode. element [4] ; 
coC4]  =  scale_factor; 

FUNCTIONS  USED 

ddphi . curvature  of  beam  at  location  x 

REFERENCES 

Harris,  C.M.  and  Crede,  C.V.  (eds).  Shock  and  Vibration  Handbook, 
2nd  edition,  McGraw-Hill,  New  York,  pp.  7-13  through  7-15,  1976. 

Thomson,  U.T.  "Theory  of  Vibration  with  Applications,"  Prentice- 
Hall,  Inc.  Englewood  Cliffs,  NJ,  1972. 
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*********************************************************************  * i 


double  *co,  x,  Z; 

{ 

double  ddphiO; 

/*  compute  and  return  the  beam  strain  at  location  x,  z.  */ 
return(  z*ddphi(  x,  co  )  ); 


double  strain2(  x,  z,  co  ) 

j *  ********************************************************************* 

NAME 


strain2 . square  of  beam  strain 

SYNOPSIS 


y  =  strain2(  x,  z,  co  );  compute  square  of  beam  strain 


double  *co, 

x. 

Y. 

*; 

DESCRIPTION 


coefficients  used  in  strain  term 
x  coordinate  of  beam 
square  of  beam  strain 
z  coordinate  of  beam 


Compute  and  square  strain  at  beam  coordinates  given  at  x  and  z. 
The  strain  is  computed  from  the  following  formula 

strain(x.z)  =  z*phi"(x) 

cotO]  =  kappa; 
coll]  =  mode. element  [2]  ; 
co£2J  =  mode. element [3]; 
cot3]  =  mode. element [4] ; 
co[4]  =  scale_fa<.tor; 
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FUNCTIONS  USED 


strain . strain  of  beam  at  location  x,  2 
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double  *co,  x,  z; 
double  strainO; 

/*  compute  and  return  the  squared  beam  strain  at  location  x,  z.  */ 
returnf  sqr(  strainf  x,  z,  co  )  )  ); 
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INTERNAL  DAMPING  TEST  RESULTS 


Results  of  the  internal  damping  test  are  presented,  together  with  a 
discussion  of  damping  measurement  enors.  Two  types  of  internal  damping 
were  observed:  amplitude-  and  frequency-dependent  damping. 


DAMPING  ERRORS  AND  CORRECTIONS 

Two  primary  sources  of  error  were  uncovered  while  collecting  the 
symmetric  cantilever  beam  damping  data.  The  first  error  involves  the 
range  of  data  collected.  Low  levels  of  stiain  and  acceleration  are 
imposed  by  the  limits  of  the  Liicke-Granato  model  (Ref  34,  44),  requiring 
sensitive  transducers  to  accurately  depict  the  low  input  levels.  In 
addition,  measurements  are  taken  in  the  first  three  flexural  modes  of 
the  6-  and  8- inch  beam  specimens;  thus,  imposing  a  large  frequency  span 
requirement  on  accelerometer  performance.  Table  E-l  lists  the  range  of 
input  data  collected  during  the  damping  tests. 


Table  E-l.  Range,  of  Accelerations  and  Strains 
for  Cantilever  Beam  Damping  Tests 


Parameter 

Range 

Strain 

10  to  100  (microstrain) 

Acceleration 

1  to  300  (in./sec^) 

Frequency 

60  to  2000  (Hz) 

E-l 


The  second  area  of  error  lies  in  the  wires  connecting  the  strain 
gages  to  their  amplifiers.  These  wires,  vibrating  with  the  beam  speci¬ 
men,  burden  the  system  with  external  damping  that  must  be  quantified  in 
order  to  extract  the  internal  damping  values  from  the  total  damping 
quantity.  Quantification  of  the  external  damping  due  to  strain  gage 
wiring  is  accomplished  by  removing  one  set  of  wires  at  a  time  and  com¬ 
puting  the  resulting  change  in  damping.  Figure  E-l  shows  the  strain 
gage  arrangement  and  indicates  the  order  in  which  the  wires  are  removed. 
Designating  the  damping  states  with  zero,  one,  two,  three  and  four  sets 

of  wires  removed  as  A^,  A ^ ,  A^,  A^,  and  a\  respectively,  indicates  a 

4  4 

damping  state,  A  ,  for  which  no  measurements  may  be  taken.  A  is 

3 

estimated  from  A  by  subtracting  the  change  in  damping  when  the  third 

3  2  3 

set  of  wires  is  removed,  that  is,  A  -  (A  -  A  )  or: 


(E-l) 


This  estimate  assumes  that  the  removal  of  the  fourth  set  of  wires  affects 
the  damping  exactly  as  the  removal  of  the  third  set  of  wires.  This  assump¬ 
tion  is  reasonable  in  light  of  the  symmetric  location  of  the  third  and 
fourth  strain  gages. 

To  determine  the  external  damping  decrement,  the  frequency-dependent 

4 

decrement,  A^,  calculated  from  Equation  37,  is  subtracted  from  A  of 
Equation  E-l.  The  external  damping  decrement  is  determined  for  a  beam 
specimen  at  the  three  fundamental  modes  and  an  average  is  taken.  The 
expression  for  the  external  damping: 


(E-2) 


where  A  and  A^  are  the  results  of  Equation  E-l  and  Equation  37  in  the 
main  text  of  this  report  for  the  i*'^1  mode,  respectively.  Figure  26  shows 
the  effect  of  Equation  E-2  on  the  frequency-dependent  damping  decrement 
for  the  first  three  modes  of  two  symmetric  cantilever  beams,  one  6  inches 
and  the  other  8  inches  long.  Figure  26  also  shows  Zener's  Theory  super- 
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imposed  on  the  experimental  data  to  demonstrate  the  error  in  using  Equation 
E-2.  The  maximum  errors  for  the  first,  second,  and  third  modes  occurred 
for  6-inch  beam  specimens  and  were  9,  36,  and  100  percent,  respectively. 

In  each  case,  the  damping  ratio  error  was  on  the  order  of  6  x  10  ^ . 

STRAIN- AMPLITUDE  DEPENDENT  DAMPING  DECREMENT 

Calculation  of  the  frequency-dependent  and  external  damping 
decrements  reduces  Equation  F-22  in  Appendix  F  to  solving  for  the 

average  strain-amplitude  dependent  damping  decrement,  A  .  In  Equation 

-  2  " 

F-17,  A  is  calculated  by  the  integration  of  e  A  (e),  where  A  (e)  is 
rl  H  n 

given  by  (Ref  44): 

Ah(e)  =  C1(C2/e)2E1(C2/e)  (E-3) 

where  and  C2  are  material  constants  and  F.j  represents  the  exponential 

integral  function  of  the  first  order. 

Figure  14  in  the  main  text  of  this  report  presents  damping  and 

strain  data  collected  from  symmetric  cantilever  beam  tests.  In  Figure 

14,  the  damping  and  strain  data  are  scaled  by  estimates  of  the  material 

constants  C.  and  C„,  respectively.  These  estimates  are  based  on  numerical 
1  2  2 

integration  of  e  A^(e)  with  A^C e )  from  Equation  E-3  and  strain  expressed 
in  the  form  of  Equation  F-24  in  Appendix  F.  A  least  squares  fit  was 
used  to  obtain  the  best  fit  of  the  data  to  the  table  generated  by  numeri¬ 
cal  integration.  An  algorithm  based  on  Neville's  Iterated  Interpolation 
was  used  to  obtain  damping  values  lying  between  tabulated  quantities. 
Converting  all  strain  measurements  to  maximum  beam  strain  resulted  in  a 

minimum  measurable  damping  ratio  of  6  x  10  ^  in  a  strain  range  of  10  to 

-3 

20  microstrain.  Estimates  of  material  constants  are  C^  =  54  x  10  and 
C2  =  144  x  10~6. 
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Figure  E-l.  Cantilever  beam  strain  gate  locations. 
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DERIVATION  OF  SYMMETRIC  CANTILEVER  BEAM  MODEL 


INTRODUCTION 

The  purpose  of  this  appendix  is  to  develop  the  equations  and  terms 
required  to  analyze  the  internal  damping  test  results.  An  equation  of 
motion  with  its  mass  and  stiffness  terms  is  presented.  The  damping 
ratio  emerges  in  a  form  directly  applicable  to  the  internal  damping  test 
approach  and  data  reduction.  Specific  damping  capacity  is  defined  and 
the  experimental  results  in  terms  of  amplitude-  and  frequency-dependent 
damping  are  expressed  together  with  external  sources  of  damping. 


EQUATION  OF  MOTION 

The  differential  equation  of  motion  for  the  symmetric  cantilever 
beam  model  is  (Ref  72): 


[EIw"(x,t)]"  +  c(x)w(x,t)  +  m(x)[w^(t)  +  w(x,t)]  =  0 


(F- 1 ) 


where  w  is  the  beam  displacement  in  the  z  direction  and  I  is  the  beam 
moment  of  inertia.  F.quation  F-l  is  rearranged  to  distinguish  the  forcing 
function : 


[EIw"(x,t)j"  +  c(x)w(x,t)  +  m(x)w(x,t)  =  -mlxlw^Ct) 


(F-2) 


F-l 


Thus,  the  forcing  function  is  expressed  as  an  inertia]  force  per  unit 
length  -m(x)w^(t).  Assuming  a  solution  in  the  form: 

w(x,t)  =  I  B.  qi(t)  <^(x)  (F-3) 

i 


where  the  second  derivative  of  the  character ist ic  shape  function  with 
respect  to  x,  $  "(x),  describes  the  curvature  in  the  vibrating  beam. 
Combining  Equation  F-7  with  the  bending  moment-bending  strain  relation¬ 
ship  (Ref  77)  produces  an  expression  for  the  amplitude-scaling  factor  in 
terms  of  a  measured  strain: 


F-2 


(  F  — 8 ) 


B.  .  -V(^L_ 
‘  h‘i  ♦/'<*„> 


where  e(x  )  is  a  measured  strain  at  location  x  and  k.  is  the  eigenvalue 
tR  o  1 

of  the  i  mode. 


DAMPINC  RATIO 


The  acceleration  expression,  on  the  right-hand  side  of  Equation 
F-4,  may  be  replaced  by  a  solitary  parameter,  z^,(t),  where: 


Zj(t)  =  -wb(t) 


&  I'  #t(x) 

Mi  J0  1 


(F-9) 


Assuming  sinusoidal  motion  with  w^(t)  =  W^sin  wt  for  the  vibrating  beam, 
it  follows  z^(t)  =  Z^sin  wt  and  z ,  =  -z^  w2  from  which  Equation  F-4  may 
be  rewritten: 


■  2  2 
q.  +  2C.w.q.  +  w.  q.  =  Z.  w  sin  wt 

l  i  i  i  i  i  l 


(  F  - 1 0  ) 


The  steady-state  solution  q.  =  sin(wt  -  0)  is  then  available  from 
inspection  to  be: 


Z . (w/w  .  ) 

g  —  . .  .  _T _  *  _  _ 

'  {[ 1  -  (w/w.)2]2  +  [2C.(w/w.)12]1/2 


(F-ll) 


where  w^  is  the  resonant  frequency  of  the  beam  and  Z  ^  is  defined  by: 


Z.  -  wb(t) 


,(x)  f’  * 
Mi  Jo  1 


(x)  dx 


When  the  beam  is  driven  at  its  resonance  frequency,  that  is,  w  =  w  , 
Equation  F-ll  reduces  to: 


CF-12) 


Rearranging  Equation  F-12  to  solve  for  the  damping  ratio,  results 

in: 


«i 


(  F  - 1 3 ) 


Determination  of  the  damping  ratio,  for  Equation  F-2  is  accomplished 

by  vibrating  the  beam  at  its  resonant  frequency,  uk,  while  measuring  the 
base  acceleration,  w^,  and  the  strain  at  some  known  location  on  the  beam, 
e(xq).  Equation  F-8  uses  z(xq)  to  solve  for  and  Equation  F-9  requires 
w^  and  uk  to  solve  for  Z^.  Using  Equation  F-13  with  and  Z^  will 
determine  the  damping  ratio. 


SPECIFIC  DAMPING  CAPACITY 


For  a  material,  the  ratio  of  energy  dissipated  per  cycle,  W^,  to 
the  peak  potential  energy  per  cycle,  U,  is  defined  as  the  specific 
damping  capacity,  A, 


A 


(F-14) 


As  developed  by  Shimanuki  and  Doi  (Ref  43),  the  average  specific  damping 
capacity  of  a  material,  subject  to  internal  friction,  is  expressed  using 
a  volume  element: 


A 


}  Aw  dV 
2  H  w  dV 


( F  - 1 5 ) 


F-4 


where  Aw  and  w  are  the  energy  dissipated  and  the  peak  potential  energy 
per  cycle  in  a  volume  element,  respectively.  Shimanuki  and  Doi  assume 
that  Aw  is  a  function  of  dislocation  damping  in  the  form  of  A  ,  the 
strain  amplitude-dependent  decrement  and  A^,  the  frequency-dependent 
decrement.  Substituting  for  Aw  and  w  in  Equation  F-15: 


A  =  Ah  +  Aj 


where 


A 


H 


I  e\(e)  dV 
f  e2  dV 


(F- 16) 


(F- 17) 


and  w  =  £  . 

In  analyzing  the  specific  damping  capacity  for  the  dynamic  model 
described  in  Equation  F-2  an  assumption  must  be  made  about  c(x).  As 
presented  in  Equation  F-16,  the  average  specific  damping  capacity  con¬ 
sists  of  two  terms:  A  which  is  strain-amplitude  dependent,  and  A 
which  is  frequency  dependent.  Thus,  neither  term  is  proportional  to 
velocity  and  c(x)  must  be  treated  as  an  equivalent  expression  of  viscous 

damping,  c  (Ref  72).  From  Equation  F-2, 
eq 


W 


d 


I 


c  w  dw 
eq 


(F- 18) 


Expanding  w  from  Equation  F-3,  Equation  F-18  becomes: 


W, ,  =  c  (  *.2(x)  dx  L  q.2(t)  dt 

d]  eq  JQ  i  f  1 1 


From  Equation  F-5: 


C  . 
eq  i 


c  (  $,2(x)dx 

eq  J0  1 


F-5 


with  the  resul  that: 


W,.  =  C  .4  q.2(t)  dt 

di  eqi  J  1 


Assuming  steady-state  harmonic  displacement  and  velocity, 


W, .  =  B  .  uwC 
d 1  l  eq l 


(F-19) 


From  Thompson  (Ref  72),  wC  .  =  2C,K,,  which  when  substituted  into 

eqi  l  i 

Equation  F-19  produces  the  final  expression  for  W^: 


W,.  =  2B.  irC.K. 

di  ill 


The  peak  potential  energy  per  cycle  in  the  i^  mode  is  given  by 


U.  =  1/2B1"K1 


Substituting  for  and  U  in  F.quation  F-14,  the  average  specific  damping 
capacity  for  the  beam  is: 


A  =  4uC . 


(F-20) 


Equating  Equations  F-16  and  F--20: 


*i[ +  4. = 


(F-21) 


Equation  F-21  assumes  no  external  damping  influences  the  damping  ratio. 
This  assumption  is  not  an  accurate  reflection  of  a  "real-world"  system 
where  attachment  points  and  connections  are  subject  to  slippage  and 
friction.  Therefore  some  external  dumping  must  be  included  in  the 


F-6 


Before  the 


system.  This  external  damping  is  represented  by  A  . 
specific  damping  capacity  can  be  determined,  A  t  must  be  estimated. 
Modifying  Equation  F-21  to  reflect  the  inclusion  of  external  damping: 

A„  +  A  =  4ir£  .  -  A  (F-22) 

HI  1  ext 

Characteristic  shape  functions  for  a  cantilever  beam  are  well- 
defined,  with  authors  using  various  formats  for  identical  expressions 
(Ref  72,  78).  This  paper  will  borrow  the  format  of  Harris  and  Crede 
(Ref  78)  to  symbolize  the  shape  function  for  a  cantilever  beam: 

<^(x)  =  (cos  k^x  -  cosh  k.x)  +  R^(sin  k^x  -  sinh  k^x)  (F-23) 

where  R^  is  a  ratio  of  nonzero  constants.  Table  F-l  displays  values  of 
k^I  and  Rj  for  the  first  three  modes  of  vibration.  Differentiating 
Equation  F-23  twice  with  respect  to  x  produces  the  beam  curvature: 

2 

0."(x)  =  -k .  [(cos  k.x  +  cosh  k.x)  +  R.(sin  k.x  +  sinh  k.x))  (F-24) 

l  l  1  l  l  i  l  i  ‘ 

Equation  F-24  is  used  in  conjunction  with  Equation  F-8  to  compute  the 
amplitude  scaling  factor,  B^,  from  experimental  data,  z(x^) . 

Table  F-l.  k,1  and  R.  for  First  Three  Modes 

l  l 

of  a  Cantilever  Beam 


Mode ,  i 

Eigenvalue,  k^I 

Ratio,  R. 

l 

1 

1.8751 

-0.7341 

2 

4.6941 

-1.0185 

3 

7.8548 

-0.9992 

F-7 
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DERIVATION  OF  FOUR-POINT  BEAM  FATIGUE  MODEL 

INTRODUCTION 

This  appendix  parallels  the  procedure  of  Appendix  F  in  the 
derivation  of  a  model  with  which  to  analyze  the  results  of  the  four- 
point  beam  fatigue  tests.  An  expression  for  the  knife  edges  supporting 
the  fatigue  specimen  is  derived  and  the  boundary  conditions  are  deter¬ 
mined.  The  equation  of  motion  and  mass  and  stiffness  terms  are  presented, 
and  the  damping  term  is  derived.  Finally,  an  empirical  expression  for 
the  damping  observed  in  4340  steel  beam  fatigue  specimens  is  provided. 


KNIFE-EDGE  SUPPORT  MODEL 

Support  for  the  four-point  beam  is  a  clamped  knife-edge  pair  (see 
Figure  10  in  the  main  text).  The  knife-edge  pair,  while  allowing  some 
rotation,  does  impose  a  moment,  M^  .  Visual  observation  of  the  knife 
edges  during  fatigue  loading  indicates  the  beam  rotates  about  the  more 
rigidly  mounted  knife  edge  of  the  pair.  This  situation  introduces  some 
axial  loading  to  the  beam  that  must  be  accounted  for  when  analyzing  the 
beam  under  fatigue  load.  However,  as  the  strains,  and  therefore  the 
rotations,  are  small  for  the  load  applied  during  damping  data  acquisition, 
it  is  assumed  the  point  of  rotation  lies  along  the  center  of  the  beam 
and  the  axial  loading  is  negligible.  The  error  associated  with  this 
assumption  is  on  the  order  of  10  ^  given  a  maximum  strain  of  100  micro 
strain  in  a  simple-supported  beam. 


G-  1 


BOUNDARY  CONDITIONS 


The  moment  imposed  by  the  knife-edge  support  is  computed  from: 
Mke=kkew’(°)  (G-l) 

where  w'(0)  is  the  beam  rotation  at  the  knife-edge  support.  As: 

Mke  =  Elw"(0) 

where  w"(0)  is  the  beam  curvature  at  the  knife-edge  support,  Equation 
G-l  provides  a  boundary  condition  in  the  form: 

Elw"(0)  =  k.  w'(0) 
ke 

The  moment  imposed  by  a  crack  at  the  beam  midpoint  is  computed  from 
(Ref  79): 

M  =  k  (a)  w'(l/2)  (G-2) 

c  c 

where  is  the  imposed  moment  and  w'(l/2)  is  the  beam  rotation  at  its 
midpoint.  The  moment  is  expressed  as  a  function  of  the  curvature  at  the 
beam  midpoint,  w"(l/2), 

Mc  =  EIw"( 1/2)  (G-3) 

Equation  G-2  and  G-3  provide  n  second  boundary  condition: 

EIw"( 1/2)  =  kc(a)  w' ( 1/2) 

The  two  remaining  boundary  conditions  are  based  on  a  beam  with  a  simple 
support  at  one  endpoint  and  guided  at  the  other.  These  boundary  condi¬ 
tions  are  given  by: 


G-2 


w(0)  =  0 


and 


w' "( 1/2)  =  0 

where  w(0)  is  the  displacement  at  the  knife-edge  support  and  w"'(l/2)  is 
the  third  derivative  of  the  beam  displacement  with  respect  to  x  at  the 
beam  midpoint. 


EQUATION  OF  MOTION 

The  equation  of  motion  for  the  four-point  beam  is  given  by: 

[EIw"(x,t)]"  +  c(x)w(x,t)  +  m(x)w(x,t)  +  m  w(l  ,t)  =  F  (1  ,t)  (G-4) 

ex  ex  ex  ex 

Employment  of  Equation  F-3  in  transforming  Equation  G-4  produces  an 
equation  of  motion  in  generalized  coordinates: 


q .  +  2( . w .q ,  + 
l  i  li 


qi 


(x) 


F 

ex 


(1 


t) 


dx 


(  G  — 5  ) 


Equations  F-5  and  F-6  apply  to  Equation  G-5;  however,  the  generalized 
mass  expression  has  been  altered  due  to  the  mass  of  the  exciter.  The 
generalized  mass  is  expressed  as: 


M. 


/: 


$  2(x)m(x)  dx  +  $,2(1  )m 
i  i  ex  ex 


The  unit  impulse  or  delta  function,  6(x-£)  (Ref  72),  models  a  point 
load  at  the  location  £  =  1  .  Assuming  F  (1  , t )  is  harmonic,  the  forcing 
function  is  separated  into  two  components,  one  x  and  the  other  t  dependent. 


G-3 


With  this  assumption,  the  forcing  function  becomes  f  6(x- lex)sin  wt, 
where  f  is  the  amplitude  scaling  factor.  Using  this  expression, 
integration  in  Equation  G-5  becomes  trivial,  with  the  result: 


q.  +  2£.W,q.  +  W  .  q. 
1  ill  11 


f 

ex  .  , .  . 

s —  ML,,)  sin  wt 

M .  i  ex 


(G-6) 


SPECIFIC  DAMPING  CAPACITY 

Treating  the  right-hand  side  of  Equation  G-6  as  the  parameter  z  (t), 
that  is, 


zi(t)  '  ST  W  sin  ut 


(G-7) 


the  relationship  between  displacement  and  specific  damping  capacity  is 

established  as  demonstrated  in  Appendix  F.  From  Equation  G-7,  it  follows 

2 

z^,(t)  =  Z^sin  wt  where  z^  =  -z  w  and: 


Z. 

l 


f 

_ex 

M 


^.(1  ) 

l  ex 


(G-8) 


Substitution  of  Z^  into  Equation  G-6  produces  Equation  F-10  of  Appendix 
F.  Thus,  the  four-point  beam  specimen  will  have  its  specific  damping 
capacity  determined  in  a  manner  similar  to  that  applied  to  the  symmetric 
cantilever  beam  damping  specimen. 

Results  from  the  4340  steel  beam  fatigue  tests  do  not  conform  to 
the  dislocation  damping  based  model.  This  appendix  develops  a  second 
damping  mechanism,  based  on  crack  surface  rotation,  that  adequately 
models  and  explains  the  observed  damping  phenomena.  This  new  mechanism 
produces  an  empirical  expression  that  is  used  in  the  four-point  beam 
fatigue  model. 
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EMPIRICAL  DAMPING  EXPRESSION 


With  the  results  of  the  4340  steel  beam  fatigue  test  suggesting  a 
strong  relationship  between  the  crack  length  and  damping,  the  decision 
was  made  to  review  previous  work  to  reconcile  the  proposed  constitutive 
relationship  with  experimental  data.  Previous  work  with  1018  steel  beams 
fP.ef  73)  included  similar  fatigue  tests.  These  tests  were  performed 
under  constant  load.  The  softer  1018  steel  forms  a  large  plastic  zone 
with  the  increasing  stress  intensity  function  associated  with  a  constant 
load.  This  large  plastic  zone  does  not  allow  the  crack  surfaces  to  touch 
as  the  beam  vibrates  to  collect  damping  data. 

Figures  16  through  18  in  the  main  text  show  a  crack  growth  rate  per 
cycle  that  is  proportional  to  the  stress  intensity  factor  raised  to  the 
fourth  power.  In  equation  form,  this  relation  is  expressed  by: 


da 

dN 


a  K, 


(G-9) 


where  ac  is  a  material  dependent  constant,  N  is  the  fatigue  cycle  number, 
and  Kj  is  the  stress  intensity  factor.  Equation  G-9  is  consistent  with 
the  dislocation  theory  of  crack  growth  as  described  by  Frost,  Marsh,  and 
Pook  (Ref  60).  Rearranging  Equation  G-9  and  integrating  over  one  cycle 
produces : 
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_ da 

Kj\a) 


a 

c 


(G- 10) 


where  the  crack  growth  in  one  cycle  is  calculated  from: 


Aa .  =  a ,  -  a ,  , 

J  J  J-l 

In  Equation  G-10,  is  assumed  to  be  a  function  of  the  crack  length. 
The  expression  for  Kj  is  taken  from  Equation  31  in  the  main  text. 
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where  b^  is  a  material -dependent  constant  (see  Figure  G-l).  For  the 
4340  steel  beam  specimens,  b^  =  2.45.  This  relationship,  for  the  three 
beam  specimens  analyzed,  resembles  the  moment  of  inertia  expression  in 
its  usage  of  beam  height  (Ref  77). 
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Figure  G-l.  Crack  stiffness  versus  crack  length  ratio. 
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